diff --git a/doc/features/aerodynamics.rst b/doc/features/aerodynamics.rst index 528a4226..26245226 100644 --- a/doc/features/aerodynamics.rst +++ b/doc/features/aerodynamics.rst @@ -31,71 +31,31 @@ In run script, users should set the values for the following aircraft design par Using OpenAeroStruct ==================== Instead of the simple drag polar model, you can use `OpenAeroStruct `_ to compute the drag. -This allows you to parametrize the aircraft wing design in more details and explore the effects of wing geometry, including taper, sweep, and twist. +This enables a more detailed parameterization of the aircraft wing. OpenAeroStruct implements the vortex-lattice method (VLM) for aerodynamics and beam-based finite element method (FEM) for structures (in the case of the aerostructural drag polar). -For more detail, please check the `documentation `_. +For more details, please check the `documentation `_. -The wing is currently limited to simple planform geometries. -Additionally, the wing does not include a tail to trim it because there is no OpenConcept weight position model with which to trim. +The aerodynamic-only model supports three types of mesh generation for planform and geometry flexibility. +The aerostructural model is currently limited to simple planform geometries. +Additionally, the wing does not include a horizontal tail to trim it. OpenConcept uses a surrogate model trained by OpenAeroStruct analyses to reduce the computational cost. -The data generation and surrogate training is automated; specifying the training grid manually may improve accuracy or decrease computational cost. +The data generation and surrogate training is automated; specifying the training grid manually may improve accuracy and decrease computational cost. VLM-based aerodynamic model: ``VLMDragPolar`` ------------------------------------------------ This model uses the vortex-lattice method (VLM) to compute the drag. The inputs to this model are the flight conditions (Mach number, altitude, dynamic pressure, lift coefficient) and aircraft design parameters. -Users should set the following design parameters and options in the run script. +The aerodynamic mesh can be defined in one of three ways: -.. list-table:: Aircraft design variables for VLM-based model - :header-rows: 1 +1. Use simple planform variables to define a trapezoidal planform. These planform variables are wing area, aspect ratio, taper ratio, and quarter chord sweep angle. - * - Variable name - - Property - - Type - * - ac|geom|wing|S_ref - - Reference wing area - - float - * - ac|geom|wing|AR - - Wing aspect ratio - - float - * - ac|geom|wing|taper - - Taper ratio - - float - * - ac|geom|wing|c4sweep - - Sweep angle at quarter chord - - float - * - ac|geom|wing|twist - - Spanwise distribution of twist, from wint tip to root. - - 1D ndarray, lendth ``num_twist`` - * - ac|aero|CD_nonwing - - Drag coefficient of components other than the wing; e.g. fuselage, - tail, interference drag, etc. - - float - * - fltcond|TempIncrement - - Temperature increment for non-standard day - - float - -.. list-table:: Options for VLM-based model - :widths: 30 50 20 - :header-rows: 1 +2. Define multiple spanwise wing sections between which a mesh is linearly interpolated. Each section is defined by its streamwise offset, chord length, and spanwise position. The whole planform is then scaled uniformly to match the desired wing area. - * - Variable name - - Property - - Type - * - ``num_x`` - - VLM mesh size (number of vertices) in chordwise direction. - - int - * - ``num_y`` - - VLM mesh size (number of vertices) in spanwise direction. - - int - * - ``num_twist`` - - Number of spanwise control points for twist distribution. - - int - -There are other advanced options, e.g., the surrogate training points in Mach-alpha-altitude space. -The default settings should work fine for these advanced options, but if you want to make changes, please refer to the source docs. +3. Directly provide an OpenAeroStruct-compatible mesh. + +More details on the inputs, outputs, and options are available in the source code documentation. Aerostructural model: ``AeroStructDragPolar`` ----------------------------------------------------- @@ -133,12 +93,12 @@ OpenConcept uses surrogate models based on OpenAeroStruct analyses to reduce the The surrogate models are trained in the 3D input space of Mach number, angle of attack, and altitude. The outputs of the surrogate models are CL and CD (and failure for ``AeroStructDragPolar``). -For more details about the surrogate models, see our `paper `_. +For more details about the surrogate models, see our `paper `_. Other models ============ -The aerodynamics module also includes a couple components that may be useful to be aware of: +The aerodynamics module also includes a couple components that may be useful: - ``StallSpeed``, which uses :math:`C_{L, \text{max}}`, aircraft weight, and wing area to compute the stall speed - ``Lift``, which computes lift force using lift coefficient, wing area, and dynamic pressure diff --git a/openconcept/__init__.py b/openconcept/__init__.py index 7863915f..6849410a 100644 --- a/openconcept/__init__.py +++ b/openconcept/__init__.py @@ -1 +1 @@ -__version__ = "1.0.2" +__version__ = "1.1.0" diff --git a/openconcept/aerodynamics/openaerostruct/__init__.py b/openconcept/aerodynamics/openaerostruct/__init__.py index b6e3e0fa..8b82ec78 100644 --- a/openconcept/aerodynamics/openaerostruct/__init__.py +++ b/openconcept/aerodynamics/openaerostruct/__init__.py @@ -1,2 +1,3 @@ -from .drag_polar import VLMDragPolar, VLMDataGen, VLM, PlanformMesh +from .mesh_gen import TrapezoidalPlanformMesh, SectionPlanformMesh, ThicknessChordRatioInterp, SectionLinearInterp +from .drag_polar import VLMDragPolar, VLMDataGen, VLM from .aerostructural import AerostructDragPolar, OASDataGen, Aerostruct, AerostructDragPolarExact diff --git a/openconcept/aerodynamics/openaerostruct/aerostructural.py b/openconcept/aerodynamics/openaerostruct/aerostructural.py index 8041383c..34c5ed84 100644 --- a/openconcept/aerodynamics/openaerostruct/aerostructural.py +++ b/openconcept/aerodynamics/openaerostruct/aerostructural.py @@ -25,7 +25,7 @@ from openaerostruct.integration.aerostruct_groups import AerostructPoint from openaerostruct.structures.spatial_beam_setup import SpatialBeamSetup from openaerostruct.structures.wingbox_group import WingboxGroup - from openconcept.aerodynamics.openaerostruct.drag_polar import PlanformMesh + from openconcept.aerodynamics.openaerostruct.mesh_gen import TrapezoidalPlanformMesh except ImportError: raise ImportError("OpenAeroStruct must be installed to use the AerostructDragPolar component") @@ -119,9 +119,9 @@ class AerostructDragPolar(om.Group): num_nodes : int Number of analysis points per mission segment (scalar, dimensionless) num_x : int - Number of points in x (streamwise) direction (scalar, dimensionless) + Number of panels in x (streamwise) direction (scalar, dimensionless) num_y : int - Number of points in y (spanwise) direction for one wing because + Number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) num_twist : int Number of spline control points for twist (scalar, dimensionless) @@ -152,8 +152,8 @@ def __init__(self, **kwargs): def initialize(self): self.options.declare("num_nodes", default=1, desc="Number of analysis points to run") - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare("num_y", default=6, desc="Number of spanwise (half wing) mesh panels") self.options.declare("num_twist", default=4, desc="Number of twist spline control points") self.options.declare("num_toverc", default=4, desc="Number of thickness to chord ratio spline control points") self.options.declare("num_skin", default=4, desc="Number of skin thickness spline control points") @@ -307,9 +307,9 @@ class OASDataGen(om.ExplicitComponent): Options ------- num_x : int - Number of points in x (streamwise) direction (scalar, dimensionless) + Number of panels in x (streamwise) direction (scalar, dimensionless) num_y : int - Number of points in y (spanwise) direction for one wing because + Number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) num_twist : int Number of spline control points for twist (scalar, dimensionless) @@ -341,8 +341,8 @@ def __init__(self, **kwargs): self.cite = CITATION def initialize(self): - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare("num_y", default=6, desc="Number of spanwise (half wing) mesh panels") self.options.declare("num_twist", default=4, desc="Number of twist spline control points") self.options.declare("num_toverc", default=4, desc="Number of thickness to chord ratio spline control points") self.options.declare("num_skin", default=4, desc="Number of skin thickness spline control points") @@ -558,9 +558,9 @@ def compute_partials(self, inputs, partials): List of spar thicknesses at control points of spline (vector, m) NOTE: length of vector is num_spar (set in options of OASDataGen) num_x: int - number of points in x (streamwise) direction (scalar, dimensionless) + number of panels in x (streamwise) direction (scalar, dimensionless) num_y: int - number of points in y (spanwise) direction for one wing because + number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) surf_dict : dict Dictionary of OpenAeroStruct surface options; any options provided here @@ -812,9 +812,9 @@ class Aerostruct(om.Group): Options ------- num_x : int - Number of points in x (streamwise) direction (scalar, dimensionless) + Number of panels in x (streamwise) direction (scalar, dimensionless) num_y : int - Number of points in y (spanwise) direction for one wing because + Number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) num_twist : int Number of spline control points for twist (scalar, dimensionless) @@ -838,8 +838,8 @@ def __init__(self, **kwargs): self.cite = CITATION def initialize(self): - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare("num_y", default=6, desc="Number of spanwise (half wing) mesh panels") self.options.declare("num_twist", default=4, desc="Number of twist spline control points") self.options.declare("num_toverc", default=4, desc="Number of thickness to chord ratio spline control points") self.options.declare("num_skin", default=4, desc="Number of skin thickness spline control points") @@ -847,8 +847,10 @@ def initialize(self): self.options.declare("surf_options", default=None, desc="Dictionary of OpenAeroStruct surface options") def setup(self): - nx = int(self.options["num_x"]) - ny = int(self.options["num_y"]) + # Number of coordinates is one more than the number of panels + nx = int(self.options["num_x"]) + 1 + ny = int(self.options["num_y"]) + 1 + n_twist = int(self.options["num_twist"]) n_skin = int(self.options["num_skin"]) n_spar = int(self.options["num_spar"]) @@ -1202,7 +1204,11 @@ def setup(self): comp.add_spline(y_cp_name="t_over_c_cp", y_interp_name="t_over_c") # Wing mesh generator - wing_group.add_subsystem("mesh_gen", PlanformMesh(num_x=nx, num_y=ny), promotes_inputs=["*"]) + wing_group.add_subsystem( + "mesh_gen", + TrapezoidalPlanformMesh(num_x=self.options["num_x"], num_y=self.options["num_y"]), + promotes_inputs=["*"], + ) # Apply twist spline to mesh wing_group.add_subsystem( @@ -1380,9 +1386,9 @@ class AerostructDragPolarExact(om.Group): num_nodes : int Number of analysis points per mission segment (scalar, dimensionless) num_x : int - Number of points in x (streamwise) direction (scalar, dimensionless) + Number of panels in x (streamwise) direction (scalar, dimensionless) num_y : int - Number of points in y (spanwise) direction for one wing because + Number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) num_twist : int Number of spline control points for twist (scalar, dimensionless) @@ -1407,8 +1413,8 @@ def __init__(self, **kwargs): def initialize(self): self.options.declare("num_nodes", default=1, desc="Number of analysis points to run") - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare("num_y", default=6, desc="Number of spanwise (half wing) mesh panels") self.options.declare("num_twist", default=4, desc="Number of twist spline control points") self.options.declare("num_toverc", default=4, desc="Number of thickness to chord ratio spline control points") self.options.declare("num_skin", default=4, desc="Number of skin thickness spline control points") @@ -1502,8 +1508,8 @@ def setup(self): def example_usage(): # Define parameters nn = 1 - num_x = 3 - num_y = 5 + num_x = 2 + num_y = 4 S = 427.8 # m^2 AR = 9.82 taper = 0.149 diff --git a/openconcept/aerodynamics/openaerostruct/drag_polar.py b/openconcept/aerodynamics/openaerostruct/drag_polar.py index adf64b81..8a17abe9 100644 --- a/openconcept/aerodynamics/openaerostruct/drag_polar.py +++ b/openconcept/aerodynamics/openaerostruct/drag_polar.py @@ -6,6 +6,12 @@ # Atmospheric calculations from openconcept.atmospherics import TemperatureComp, PressureComp, DensityComp, SpeedOfSoundComp +from openconcept.aerodynamics.openaerostruct import ( + TrapezoidalPlanformMesh, + SectionPlanformMesh, + ThicknessChordRatioInterp, + SectionLinearInterp, +) # Progress bar progress_bar = True @@ -39,15 +45,26 @@ class VLMDragPolar(om.Group): """ Drag polar generated using OpenAeroStruct's vortex lattice method and a surrogate - model to decrease the computational cost. + model to decrease the computational cost. It allows for three planform definitions. + The first is a trapezoidal planform defined by aspect ratio, taper, and sweep. The + second builds a planform from sections definitions, each defined by a streamwise + offset of the leading edge, spanwise position, and chord length. The mesh is linearly + lofted between the sections. Finally, this component can take in a mesh directly to + enable more flexibility. This component enables twisting of the mesh as well. Notes ----- Twist is ordered starting at the tip and moving to the root; a twist - of [-1, 0, 1] would have a tip twist of -1 deg and root twist of 1 deg + of [-1, 0, 1] would have a tip twist of -1 deg and root twist of 1 deg. + This is the same ordering as for the section definitions if the geometry + option is set to \"section\". Set the OMP_NUM_THREADS environment variable to 1 for much better parallel training performance! + Make sure that the wing reference area provided to this group is consistent both with the reference + area used by OpenAeroStruct (projected area by default) AND the reference area used by OpenConcept + to compute the lift coefficient it provides. + Inputs ------ fltcond|CL : float @@ -60,15 +77,6 @@ class VLMDragPolar(om.Group): Dynamic pressure (vector, Pascals) ac|geom|wing|S_ref : float Full planform area (scalar, m^2) - ac|geom|wing|AR : float - Aspect ratio (scalar, dimensionless) - ac|geom|wing|taper : float - Taper ratio (must be >0 and <=1); tip chord / root chord (scalar, dimensionless) - ac|geom|wing|c4sweep : float - Quarter chord sweep angle (scalar, degrees) - ac|geom|wing|twist : float - List of twist angles at control points of spline (vector, degrees) - NOTE: length of vector is num_twist (set in options), NOT num_nodes ac|aero|CD_nonwing : float Drag coefficient of components other than the wing; e.g. fuselage, tail, interference drag, etc.; this value is simply added to the @@ -76,27 +84,70 @@ class VLMDragPolar(om.Group): fltcond|TempIncrement : float Temperature increment for non-standard day (scalar, degC) NOTE: fltcond|TempIncrement is a scalar in this component but a vector in OC. \ - This will be the case for the forseeable future because of the way the \ - VLMDataGen component is set up. To make it work, TempIncrement would \ - need to be an input to the surrogate, which is not worth the extra \ - training cost (at minimum a 2x increase). + This is the case because of the way the VLMDataGen component is set up. To \ + make it work, TempIncrement would need to be an input to the surrogate, \ + which is not worth the extra training cost (at minimum a 2x increase). + + If geometry option is \"trapezoidal\" (the default) + - **ac|geom|wing|AR** *(float)* - Aspect ratio (scalar, dimensionless) + - **ac|geom|wing|taper** *(float)* - Taper ratio (must be >0 and <=1); tip chord / root chord (scalar, dimensionless) + - **ac|geom|wing|c4sweep** *(float)* - Quarter chord sweep angle (scalar, degrees) + - **ac|geom|wing|toverc** *(float)* - Wing tip and wing root airfoil thickness to chord ratio in that order; panel + thickness to chord ratios are linearly interpolated between the tip and root (vector of length 2, dimensionless) + - **ac|geom|wing|twist** *(float)* - Twist at spline control points, ordered from tip to root (vector of length num_twist, degrees) + + If geometry option is \"section\" (despite inputs in m, they are scaled to provided planform area) + - **ac|geom|wing|x_LE_sec** *(float)* - Streamwise offset of the section's leading edge, starting with the outboard + section (wing tip) and moving inboard toward the root (vector of length + num_sections, m) + - **ac|geom|wing|y_sec** *(float)* - Spanwise location of each section, starting with the outboard section (wing + tip) at the MOST NEGATIVE y value and moving inboard (increasing y value) + toward the root; the user does not provide a value for the root because it + is always 0.0 (vector of length num_sections - 1, m) + - **ac|geom|wing|chord_sec** *(float)* - Chord of each section, starting with the outboard section (wing tip) and + moving inboard toward the root (vector of length num_sections, m) + - **ac|geom|wing|toverc** *(float)* - Thickness to chord ratio of the airfoil at each defined section starting at the wing tip + and moving wing root airfoil; panel thickness to chord ratios is linearly interpolated (vector of length num_sections, m) + - **ac|geom|wing|twist** *(float)* - Twist at each section, ordered from tip to root (vector of length num_sections, degrees) + + If geometry option is \"mesh\" + - **ac|geom|wing|OAS_mesh** *(float)* - OpenAeroStruct 3D mesh (num_x + 1 x num_y + 1 x 3 vector, m) + - **ac|geom|wing|toverc** *(float)* - Thickness to chord ratio of each streamwise strip of panels ordered from + wing tip to wing root (vector of length num_y, dimensionless) Outputs ------- drag : float Drag force (vector, Newtons) + twisted_mesh : float + OpenAeroStruct mesh that has been twisted according to twist inputs; only an output if geometry options + is set to \"trapezoidal\" or \"section\" (num_x + 1 x num_y + 1 x 3 vector, m) Options ------- num_nodes : int Number of analysis points per mission segment (scalar, dimensionless) + geometry : float + Choose the geometry parameterization from the following options (by default trapezoidal): + + - "trapezoidal": Create a trapezoidal mesh based on apsect ratio, taper, and sweep + - "section": Create a mesh lofted between defined sections + - "mesh": Pass in a mesh directly to this component + num_x : int - Number of points in x (streamwise) direction (scalar, dimensionless) + Number of panels in x (streamwise) direction (scalar, dimensionless) num_y : int - Number of points in y (spanwise) direction for one wing because - uses symmetry (scalar, dimensionless) + Number of panels in y (spanwise) direction for half wing. If geometry option + is set to \"section\", this value represents the number of panels between each + pair of sections. In that case, it can also be a vector with potentially a different + number of panels between each pair of sections (scalar or vector, dimensionless) + num_sections : int + Only if geometry option is \"section\", this represents the number of spanwise sections + to define planform shape. This parameter is ignored for other geometry options (scalar, dimensionless) num_twist : int - Number of spline control points for twist (scalar, dimensionless) + Number of spline control points for twist, only used if geometry is set to \"trapezoidal\" because + \"mesh\" linearly interpolates twist between sections and \"mesh\" does not provide twist + functionality (scalar, dimensionless) alpha_train : list or ndarray List of angle of attack values at which to train the model (ndarray, degrees) Mach_train : list or ndarray @@ -107,9 +158,8 @@ class VLMDragPolar(om.Group): Dictionary of OpenAeroStruct surface options; any options provided here will override the default ones; see the OpenAeroStruct documentation for more information. Because the geometry transformations are excluded in this model (to simplify the interface), - the _cp options are not supported. The input ac|geom|wing|twist is the same - as modifying the twist_cp option in the surface dictionary. The mesh geometry modification - is limited to adjusting the input parameters to this component. + the _cp options are not supported. The t_over_c option is also removed since + it is provided instead as an input. """ def __init__(self, **kwargs): @@ -118,49 +168,159 @@ def __init__(self, **kwargs): def initialize(self): self.options.declare("num_nodes", default=1, desc="Number of analysis points to run") - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") + self.options.declare( + "geometry", + default="trapezoidal", + values=["trapezoidal", "section", "mesh"], + desc="OpenAeroStruct mesh parameterization", + ) + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare( + "num_y", default=6, types=(int, list, tuple, np.ndarray), desc="Number of spanwise (half wing) mesh panels" + ) + self.options.declare( + "num_sections", default=2, types=int, desc="Number of sections along the half span to define" + ) self.options.declare("num_twist", default=4, desc="Number of twist spline control points") self.options.declare( "Mach_train", default=np.array([0.1, 0.3, 0.45, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9]), + types=(list, np.ndarray), desc="List of Mach number training values (dimensionless)", ) self.options.declare( - "alpha_train", default=np.linspace(-10, 15, 6), desc="List of angle of attack training values (degrees)" + "alpha_train", + default=np.linspace(-10, 15, 6), + types=(list, np.ndarray), + desc="List of angle of attack training values (degrees)", ) self.options.declare( - "alt_train", default=np.linspace(0, 12e3, 4), desc="List of altitude training values (meters)" + "alt_train", + default=np.linspace(0, 12e3, 4), + types=(list, np.ndarray), + desc="List of altitude training values (meters)", + ) + self.options.declare( + "surf_options", default={}, types=(dict), desc="Dictionary of OpenAeroStruct surface options" ) - self.options.declare("surf_options", default={}, desc="Dictionary of OpenAeroStruct surface options") def setup(self): nn = self.options["num_nodes"] + nx = self.options["num_x"] + ny = self.options["num_y"] + n_twist = self.options["num_twist"] + geo = self.options["geometry"] n_alpha = self.options["alpha_train"].size n_Mach = self.options["Mach_train"].size n_alt = self.options["alt_train"].size + # Get the total number of y coordinates if it's not a scalar + if geo == "section" and isinstance(ny, int): + ny_tot = (self.options["num_sections"] - 1) * ny + 1 + elif not isinstance(ny, int): + ny_tot = np.sum(ny) + 1 + else: + ny_tot = ny + 1 + + # Generate the mesh if geometry parameterization says we should and twist it + if geo == "trapezoidal": + if not isinstance(ny, int): + raise ValueError("The num_y option must be an integer if the geometry option is trapezoidal") + + self.add_subsystem( + "gen_mesh", + TrapezoidalPlanformMesh(num_x=nx, num_y=ny), + promotes_inputs=[ + ("S", "ac|geom|wing|S_ref"), + ("AR", "ac|geom|wing|AR"), + ("taper", "ac|geom|wing|taper"), + ("sweep", "ac|geom|wing|c4sweep"), + ], + promotes_outputs=[("mesh", "ac|geom|wing|OAS_mesh")], + ) + + # Compute the twist at mesh y-coordinates + comp = self.add_subsystem( + "twist_bsp", + om.SplineComp( + method="bsplines", + x_interp_val=np.linspace(0, 1, ny_tot), + num_cp=n_twist, + interp_options={"order": min(n_twist, 4)}, + ), + promotes_inputs=[("twist_cp", "ac|geom|wing|twist")], + ) + comp.add_spline(y_cp_name="twist_cp", y_interp_name="twist", y_units="deg") + self.set_input_defaults("ac|geom|wing|twist", np.zeros(n_twist), units="deg") + + # Apply twist spline to mesh + self.add_subsystem( + "twist_mesh", + Rotate(val=np.zeros(ny_tot), mesh_shape=(nx + 1, ny_tot, 3), symmetry=True), + promotes_inputs=[("in_mesh", "ac|geom|wing|OAS_mesh")], + promotes_outputs=[("mesh", "twisted_mesh")], + ) + self.connect("twist_bsp.twist", "twist_mesh.twist") + elif geo == "section": + self.add_subsystem( + "gen_mesh", + SectionPlanformMesh(num_x=nx, num_y=ny, num_sections=self.options["num_sections"], scale_area=True), + promotes_inputs=[ + ("S", "ac|geom|wing|S_ref"), + ("x_LE_sec", "ac|geom|wing|x_LE_sec"), + ("y_sec", "ac|geom|wing|y_sec"), + ("chord_sec", "ac|geom|wing|chord_sec"), + ], + promotes_outputs=[("mesh", "ac|geom|wing|OAS_mesh")], + ) + + # Compute the twist at mesh y-coordinates + self.add_subsystem( + "twist_sec", + SectionLinearInterp(num_y=ny, num_sections=self.options["num_sections"], units="deg", cos_spacing=True), + promotes_inputs=[("property_sec", "ac|geom|wing|twist")], + ) + self.set_input_defaults("ac|geom|wing|twist", np.zeros(self.options["num_sections"]), units="deg") + + # Apply twist spline to mesh + self.add_subsystem( + "twist_mesh", + Rotate(val=np.zeros(ny_tot), mesh_shape=(nx + 1, ny_tot, 3), symmetry=True), + promotes_inputs=[("in_mesh", "ac|geom|wing|OAS_mesh")], + promotes_outputs=[("mesh", "twisted_mesh")], + ) + self.connect("twist_sec.property_node", "twist_mesh.twist") + + # Interpolate the thickness to chord ratios for each panel + if geo in ["trapezoidal", "section"]: + n_sections = self.options["num_sections"] if geo == "section" else 2 + self.add_subsystem( + "t_over_c_interp", + ThicknessChordRatioInterp(num_y=ny, num_sections=n_sections), + promotes_inputs=[("section_toverc", "ac|geom|wing|toverc")], + ) + + # Inputs to promote from the calls to OpenAeroStruct (if geometry is "mesh", + # promote the thickness-to-chord ratio directly) + VLM_promote_inputs = ["ac|aero|CD_nonwing", "fltcond|TempIncrement"] + if geo == "mesh": + VLM_promote_inputs += ["ac|geom|wing|toverc", "ac|geom|wing|OAS_mesh"] + else: + self.connect("t_over_c_interp.panel_toverc", "training_data.ac|geom|wing|toverc") + self.connect("twisted_mesh", "training_data.ac|geom|wing|OAS_mesh") + # Training data self.add_subsystem( "training_data", VLMDataGen( - num_x=self.options["num_x"], - num_y=self.options["num_y"], - num_twist=self.options["num_twist"], + num_x=nx, + num_y=ny_tot - 1, alpha_train=self.options["alpha_train"], Mach_train=self.options["Mach_train"], alt_train=self.options["alt_train"], surf_options=self.options["surf_options"], ), - promotes_inputs=[ - "ac|geom|wing|S_ref", - "ac|geom|wing|AR", - "ac|geom|wing|taper", - "ac|geom|wing|c4sweep", - "ac|geom|wing|twist", - "ac|aero|CD_nonwing", - "fltcond|TempIncrement", - ], + promotes_inputs=VLM_promote_inputs, ) # Surrogate model @@ -202,28 +362,26 @@ def setup(self): ) self.connect("aero_surrogate.CD", "drag_calc.CD") + # Set input defaults for inputs promoted from different places with different values + self.set_input_defaults("ac|geom|wing|S_ref", 1.0, units="m**2") + class VLMDataGen(om.ExplicitComponent): """ Generates a grid of OpenAeroStruct lift and drag data to train a surrogate model. The grid is defined by the options and the - planform geometry by the inputs. This component will only recalculate - the lift and drag grid when the planform shape changes. All VLMDataGen + mesh geometry by the input. This component will only recalculate + the lift and drag grid when the mesh changes. All VLMDataGen components in the model must use the same training points and surf_options. Inputs ------ - ac|geom|wing|S_ref : float - Full planform area (scalar, m^2) - ac|geom|wing|AR : float - Aspect ratio (scalar, dimensionless) - ac|geom|wing|taper : float - Taper ratio (must be >0 and <=1); tip chord / root chord (scalar, dimensionless) - ac|geom|wing|c4sweep : float - Quarter chord sweep angle (scalar, degrees) - ac|geom|wing|twist : float - List of twist angles at control points of spline (vector, degrees) - NOTE: length of vector is num_twist (set in options) + ac|geom|wing|OAS_mesh: float + OpenAeroStruct 3D mesh (num_x + 1 x num_y + 1 x 3 vector, m) + ac|geom|wing|toverc : float + Thickness to chord ratio of each streamwise strip of panels ordered from wing + tip to wing root; used for the viscous and wave drag calculations + (vector of length num_y, dimensionless) ac|aero|CD_nonwing : float Drag coefficient of components other than the wing; e.g. fuselage, tail, interference drag, etc.; this value is simply added to the @@ -241,12 +399,10 @@ class VLMDataGen(om.ExplicitComponent): Options ------- num_x : int - Number of points in x (streamwise) direction (scalar, dimensionless) + Number of panels in x (streamwise) direction (scalar, dimensionless) num_y : int - Number of points in y (spanwise) direction for one wing because + Number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) - num_twist : int - Number of spline control points for twist (scalar, dimensionless) Mach_train : list or ndarray List of Mach numbers at which to train the model (ndarray, dimensionless) alpha_train : list or ndarray @@ -257,9 +413,8 @@ class VLMDataGen(om.ExplicitComponent): Dictionary of OpenAeroStruct surface options; any options provided here will override the default ones; see the OpenAeroStruct documentation for more information. Because the geometry transformations are excluded in this model (to simplify the interface), - the _cp options are not supported. The input ac|geom|wing|twist is the same - as modifying the twist_cp option in the surface dictionary. The mesh geometry modification - is limited to adjusting the input parameters to this component. + the _cp options are not supported. The t_over_c option is also removed since + it is provided instead as an input. """ def __init__(self, **kwargs): @@ -267,9 +422,8 @@ def __init__(self, **kwargs): self.cite = CITATION def initialize(self): - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") - self.options.declare("num_twist", default=4, desc="Number of twist spline control points") + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare("num_y", default=6, desc="Number of spanwise (half wing) mesh panels") self.options.declare( "Mach_train", default=np.array([0.1, 0.3, 0.45, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9]), @@ -284,16 +438,11 @@ def initialize(self): self.options.declare("surf_options", default={}, desc="Dictionary of OpenAeroStruct surface options") def setup(self): - self.add_input("ac|geom|wing|S_ref", units="m**2") - self.add_input("ac|geom|wing|AR") - self.add_input("ac|geom|wing|taper") - self.add_input("ac|geom|wing|c4sweep", units="deg") - self.add_input( - "ac|geom|wing|twist", - val=np.zeros(self.options["num_twist"]), - shape=(self.options["num_twist"],), - units="deg", - ) + nx = self.options["num_x"] + ny = self.options["num_y"] + + self.add_input("ac|geom|wing|OAS_mesh", units="m", shape=(nx + 1, ny + 1, 3)) + self.add_input("ac|geom|wing|toverc", shape=(ny,), val=0.12) self.add_input("ac|aero|CD_nonwing", val=0.0) self.add_input("fltcond|TempIncrement", val=0.0, units="degC") @@ -327,11 +476,8 @@ def setup(self): VLMDataGen.surf_options = deepcopy(self.options["surf_options"]) # Generate grids and default cached values for training inputs and outputs - VLMDataGen.S = -np.ones(1) - VLMDataGen.AR = -np.ones(1) - VLMDataGen.taper = -np.ones(1) - VLMDataGen.c4sweep = -np.ones(1) - VLMDataGen.twist = -1 * np.ones((self.options["num_twist"],)) + VLMDataGen.mesh = -np.ones((nx + 1, ny + 1, 3)) + VLMDataGen.toverc = -np.ones(ny) VLMDataGen.temp_incr = -42 * np.ones(1) VLMDataGen.Mach, VLMDataGen.alpha, VLMDataGen.alt = np.meshgrid( self.options["Mach_train"], self.options["alpha_train"], self.options["alt_train"], indexing="ij" @@ -341,35 +487,25 @@ def setup(self): VLMDataGen.partials = None def compute(self, inputs, outputs): - S = inputs["ac|geom|wing|S_ref"] - AR = inputs["ac|geom|wing|AR"] - taper = inputs["ac|geom|wing|taper"] - sweep = inputs["ac|geom|wing|c4sweep"] - twist = inputs["ac|geom|wing|twist"] + mesh = inputs["ac|geom|wing|OAS_mesh"] + toverc = inputs["ac|geom|wing|toverc"] CD_nonwing = inputs["ac|aero|CD_nonwing"] temp_incr = inputs["fltcond|TempIncrement"] # If the inputs are unchaged, use the previously calculated values tol = 1e-13 # floating point comparison tolerance if ( - np.abs(S - VLMDataGen.S) < tol - and np.abs(AR - VLMDataGen.AR) < tol - and np.abs(taper - VLMDataGen.taper) < tol - and np.abs(sweep - VLMDataGen.c4sweep) < tol - and np.all(np.abs(twist - VLMDataGen.twist) < tol) + np.all(np.abs(mesh - VLMDataGen.mesh) < tol) + and np.all(np.abs(toverc - VLMDataGen.toverc) < tol) and np.abs(temp_incr - VLMDataGen.temp_incr) < tol ): outputs["CL_train"] = VLMDataGen.CL outputs["CD_train"] = VLMDataGen.CD + CD_nonwing return - print(f"S = {S}; AR = {AR}; taper = {taper}; sweep = {sweep}; twist = {twist}; temp_incr = {temp_incr}") # Copy new values to cached ones - VLMDataGen.S[:] = S - VLMDataGen.AR[:] = AR - VLMDataGen.taper[:] = taper - VLMDataGen.c4sweep[:] = sweep - VLMDataGen.twist[:] = twist + VLMDataGen.mesh[:, :, :] = mesh + VLMDataGen.toverc[:] = toverc VLMDataGen.temp_incr[:] = temp_incr # Compute new training values @@ -378,11 +514,8 @@ def compute(self, inputs, outputs): train_in["alpha_grid"] = VLMDataGen.alpha train_in["alt_grid"] = VLMDataGen.alt train_in["TempIncrement"] = temp_incr - train_in["S_ref"] = S - train_in["AR"] = AR - train_in["taper"] = taper - train_in["c4sweep"] = sweep - train_in["twist"] = twist + train_in["mesh"] = mesh + train_in["toverc"] = toverc train_in["num_x"] = self.options["num_x"] train_in["num_y"] = self.options["num_y"] @@ -417,29 +550,23 @@ def compute_partials(self, inputs, partials): Altitude (3D meshgrid 'ij'-style ndarray, m) TempIncrement : float Temperature increment for non-standard day (scalar, degC) - S_ref : float - Wing planform area (scalar, m^2) - AR : float - Wing aspect ratio (scalar, dimensionless) - taper : float - Wing taper ratio (scalar, dimensionless) - c4sweep : float - Wing sweep measured at quarter chord (scalar, degrees) - twist : float - List of twist angles at control points of spline (vector, degrees) - NOTE: length of vector is num_twist (set in options of VLMDataGen) + mesh: ndarray + OpenAeroStruct 3D mesh (num_x + 1 x num_y + 1 x 3 ndarray, m) + toverc : float + Thickness to chord ratio of each streamwise strip of panels ordered from wing + tip to wing root; used for the viscous and wave drag calculations + (vector of length num_y, dimensionless) num_x: int - number of points in x (streamwise) direction (scalar, dimensionless) + number of panels in x (streamwise) direction (scalar, dimensionless) num_y: int - number of points in y (spanwise) direction for one wing because + number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) surf_dict : dict Dictionary of OpenAeroStruct surface options; any options provided here will override the default ones; see the OpenAeroStruct documentation for more information. Because the geometry transformations are excluded in this model (to simplify the interface), - the _cp options are not supported. The input ac|geom|wing|twist is the same - as modifying the twist_cp option in the surface dictionary. The mesh geometry modification - is limited to adjusting the input parameters to this component. + the _cp options are not supported. The t_over_c option is also removed since + it is provided instead as an input. Returns ------- @@ -469,7 +596,7 @@ def compute_training_data(inputs, surf_dict=None): ] ).T.tolist() inputs_to_send = {"surf_dict": surf_dict} - keys = ["TempIncrement", "S_ref", "AR", "taper", "c4sweep", "twist", "num_x", "num_y"] + keys = ["TempIncrement", "mesh", "toverc", "num_x", "num_y"] for key in keys: inputs_to_send[key] = inputs[key] for row in test_points: @@ -486,31 +613,22 @@ def compute_training_data(inputs, surf_dict=None): CL = np.zeros(inputs["Mach_number_grid"].shape) CD = np.zeros(inputs["Mach_number_grid"].shape) jac_num_rows = inputs["Mach_number_grid"].size # product of array dimensions + mesh_jac_num_cols = inputs["mesh"].size + num_y_panels = inputs["mesh"].shape[1] - 1 of = ["CL_train", "CD_train"] - wrt = [ - "ac|geom|wing|S_ref", - "ac|geom|wing|AR", - "ac|geom|wing|taper", - "ac|geom|wing|c4sweep", - "ac|geom|wing|twist", - "fltcond|TempIncrement", - ] partials = {} for f in of: - for u in wrt: - if u == "ac|geom|wing|twist": - partials[f, u] = np.zeros((jac_num_rows, inputs["twist"].size)) - else: - partials[f, u] = np.zeros((jac_num_rows, 1)) + partials[f, "ac|geom|wing|OAS_mesh"] = np.zeros((jac_num_rows, mesh_jac_num_cols)) + partials[f, "ac|geom|wing|toverc"] = np.zeros((jac_num_rows, num_y_panels)) + partials[f, "fltcond|TempIncrement"] = np.zeros((jac_num_rows, 1)) data = {"CL": CL, "CD": CD, "partials": partials} # Transfer data into output data structure the proper format for i in range(len(out)): data["CL"][np.unravel_index(i, inputs["Mach_number_grid"].shape)] = out[i]["CL"] data["CD"][np.unravel_index(i, inputs["Mach_number_grid"].shape)] = out[i]["CD"] - for f in of: - for u in wrt: - data["partials"][f, u][i] = out[i]["partials"][f, u] + for key in data["partials"].keys(): + data["partials"][key][i] = out[i]["partials"][key] print(f" ...done in {time() - t_start} sec\n") @@ -530,21 +648,17 @@ def compute_aerodynamic_data(point): VLM( num_x=inputs["num_x"], num_y=inputs["num_y"], - num_twist=inputs["twist"].size, surf_options=inputs["surf_dict"], ), promotes=["*"], ) - # Set design variables - p.model.set_input_defaults("fltcond|TempIncrement", val=inputs["TempIncrement"], units="degC") - p.model.set_input_defaults("ac|geom|wing|S_ref", val=inputs["S_ref"], units="m**2") - p.model.set_input_defaults("ac|geom|wing|AR", val=inputs["AR"]) - p.model.set_input_defaults("ac|geom|wing|taper", val=inputs["taper"]) - p.model.set_input_defaults("ac|geom|wing|c4sweep", val=inputs["c4sweep"], units="deg") - p.model.set_input_defaults("ac|geom|wing|twist", val=inputs["twist"], units="deg") - p.setup() + p.setup(mode="rev") + # Set variables + p.set_val("fltcond|TempIncrement", val=inputs["TempIncrement"], units="degC") + p.set_val("ac|geom|wing|OAS_mesh", val=inputs["mesh"], units="m") + p.set_val("ac|geom|wing|toverc", val=inputs["toverc"]) p.set_val("fltcond|M", point[0]) p.set_val("fltcond|alpha", point[1], units="deg") p.set_val("fltcond|h", point[2], units="m") @@ -556,15 +670,10 @@ def compute_aerodynamic_data(point): # Compute derivatives of = ["fltcond|CL", "fltcond|CD"] of_out = ["CL_train", "CD_train"] - wrt = [ - "ac|geom|wing|S_ref", - "ac|geom|wing|AR", - "ac|geom|wing|taper", - "ac|geom|wing|c4sweep", - "ac|geom|wing|twist", - "fltcond|TempIncrement", - ] + wrt = ["ac|geom|wing|OAS_mesh", "ac|geom|wing|toverc", "fltcond|TempIncrement"] + deriv = p.compute_totals(of, wrt) + for n, f in enumerate(of): for u in wrt: output["partials"][of_out[n], u] = np.copy(deriv[f, u]) @@ -574,12 +683,12 @@ def compute_aerodynamic_data(point): class VLM(om.Group): """ - Computes lift and drag using OpenAeroStruct's vortex lattice implementation. + Computes lift, drag, and other forces using OpenAeroStruct's vortex lattice implementation. + This group contains a Reynolds number calculation that uses a linear approximation of dynamic + viscosity. It is accurate up through ~33k ft and probably ok up to 40k ft, but not much further. - Notes - ----- - Twist is ordered starting at the tip and moving to the root; a twist - of [-1, 0, 1] would have a tip twist of -1 deg and root twist of 1 deg + Thickness to chord ratio of the airfoils is taken from the input, not the value in + the surf_options dict (which is ignored). Inputs ------ @@ -591,17 +700,12 @@ class VLM(om.Group): Altitude (scalar, m) fltcond|TempIncrement : float Temperature increment for non-standard day (scalar, degC) - ac|geom|wing|S_ref : float - Wing planform area (scalar, m^2) - ac|geom|wing|AR : float - Wing aspect ratio (scalar, dimensionless) - ac|geom|wing|taper : float - Wing taper ratio (scalar, dimensionless) - ac|geom|wing|c4sweep : float - Wing sweep measured at quarter chord (scalar, degrees) - ac|geom|wing|twist : float - List of twist angles at control points of spline (vector, degrees) - NOTE: length of vector is num_twist (set in options) + ac|geom|wing|OAS_mesh: float + OpenAeroStruct 3D mesh (num_x + 1 x num_y + 1 x 3 vector, m) + ac|geom|wing|toverc : float + Thickness to chord ratio of each streamwise strip of panels ordered from wing + tip to wing root; used for the viscous and wave drag calculations + (vector of length num_y, dimensionless) Outputs ------- @@ -609,23 +713,30 @@ class VLM(om.Group): Lift coefficient of wing (scalar, dimensionless) fltcond|CD : float Drag coefficient of wing (scalar, dimensionless) + fltcond|CDi : float + Induced drag component (scalar, dimensionless) + fltcond|CDv : float + Viscous drag component (scalar, dimensionless) + fltcond|CDw : float + Wave drag component (scalar, dimensionless) + sectional_CL : float + Sectional lift coefficient for each chordwise panel strip (vector of length num_y, dimensionless) + panel_forces : float + Force from VLM for each panel (vector of size (num_x x num_y x 3), N) Options ------- num_x : int - Number of points in x (streamwise) direction (scalar, dimensionless) + Number of panels in x (streamwise) direction (scalar, dimensionless) num_y : int - Number of points in y (spanwise) direction for one wing because + Number of panels in y (spanwise) direction for one wing because uses symmetry (scalar, dimensionless) - num_twist : int - Number of spline control points for twist (scalar, dimensionless) surf_options : dict Dictionary of OpenAeroStruct surface options; any options provided here will override the default ones; see the OpenAeroStruct documentation for more information. Because the geometry transformations are excluded in this model (to simplify the interface), - the _cp options are not supported. The input ac|geom|wing|twist is the same - as modifying the twist_cp option in the surface dictionary. The mesh geometry modification - is limited to adjusting the input parameters to this component. + the _cp options are not supported. The t_over_c option is also removed since + it is provided instead as an input. """ def __init__(self, **kwargs): @@ -633,45 +744,14 @@ def __init__(self, **kwargs): self.cite = CITATION def initialize(self): - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") - self.options.declare("num_twist", default=4, desc="Number of twist spline control points") + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare("num_y", default=6, desc="Number of spanwise (half wing) mesh panels") self.options.declare("surf_options", default=None, desc="Dictionary of OpenAeroStruct surface options") def setup(self): - nx = int(self.options["num_x"]) - ny = int(self.options["num_y"]) - n_twist = int(self.options["num_twist"]) - - # ================================================================= - # Set up mesh - # ================================================================= - self.add_subsystem( - "mesh", - PlanformMesh(num_x=nx, num_y=ny), - promotes_inputs=[ - ("S", "ac|geom|wing|S_ref"), - ("AR", "ac|geom|wing|AR"), - ("taper", "ac|geom|wing|taper"), - ("sweep", "ac|geom|wing|c4sweep"), - ], - ) - - # Add bspline component for twist - x_interp = np.linspace(0.0, 1.0, ny) - comp = self.add_subsystem( - "twist_bsp", - om.SplineComp( - method="bsplines", x_interp_val=x_interp, num_cp=n_twist, interp_options={"order": min(n_twist, 4)} - ), - promotes_inputs=[("twist_cp", "ac|geom|wing|twist")], - ) - comp.add_spline(y_cp_name="twist_cp", y_interp_name="twist", y_units="deg") - - # Apply twist spline to mesh - self.add_subsystem("twist_mesh", Rotate(val=np.zeros(ny), mesh_shape=(nx, ny, 3), symmetry=True)) - self.connect("twist_bsp.twist", "twist_mesh.twist") - self.connect("mesh.mesh", "twist_mesh.in_mesh") + # Number of coordinates is one more than the number of panels + nx = int(self.options["num_x"]) + 1 + ny = int(self.options["num_y"]) + 1 # ================================================================= # Compute atmospheric and fluid properties @@ -716,7 +796,7 @@ def setup(self): np.linspace(0, 1, nx), np.linspace(-1, 0, ny), indexing="ij" ) - surf_dict = { + self.surf_dict = { "name": "wing", "mesh": dummy_mesh, # this must be defined # because the VLMGeometry component uses the shape of the mesh in this @@ -746,185 +826,45 @@ def setup(self): # Overwrite any options in the surface dict with those provided in the options if self.options["surf_options"] is not None: for key in self.options["surf_options"]: - surf_dict[key] = self.options["surf_options"][key] + self.surf_dict[key] = self.options["surf_options"][key] self.add_subsystem( "aero_point", - AeroPoint(surfaces=[surf_dict]), - promotes_inputs=[("Mach_number", "fltcond|M"), ("alpha", "fltcond|alpha")], + AeroPoint(surfaces=[self.surf_dict]), + promotes_inputs=[ + ("Mach_number", "fltcond|M"), + ("alpha", "fltcond|alpha"), + (f"{self.surf_dict['name']}_perf.t_over_c", "ac|geom|wing|toverc"), + ], promotes_outputs=[ - (f"{surf_dict['name']}_perf.CD", "fltcond|CD"), - (f"{surf_dict['name']}_perf.CL", "fltcond|CL"), + (f"{self.surf_dict['name']}_perf.CD", "fltcond|CD"), + (f"{self.surf_dict['name']}_perf.CL", "fltcond|CL"), + (f"{self.surf_dict['name']}_perf.CDi", "fltcond|CDi"), + (f"{self.surf_dict['name']}_perf.CDv", "fltcond|CDv"), + (f"{self.surf_dict['name']}_perf.CDw", "fltcond|CDw"), + (f"{self.surf_dict['name']}_perf.Cl", "sectional_CL"), + (f"aero_states.{self.surf_dict['name']}_sec_forces", "panel_forces"), ], ) - self.connect( - "twist_mesh.mesh", - [f"aero_point.{surf_dict['name']}.def_mesh", f"aero_point.aero_states.{surf_dict['name']}_def_mesh"], - ) self.connect("airspeed.Utrue", "aero_point.v") self.connect("density.fltcond|rho", "aero_point.rho") self.connect("Re_calc.re", "aero_point.re") - # Set input defaults for inputs that go to multiple locations - self.set_input_defaults("fltcond|M", 0.1) - self.set_input_defaults("fltcond|alpha", 0.0) - - # Set the thickness to chord ratio for wave and viscous drag calculation. - # It must have a thickness to chord ratio for each panel, so there must be - # ny-1 elements. Allow either one value (and duplicate it ny-1 times) or - # an array of length ny-1, but nothing else. - # NOTE: for aerostructural cases, this should be a design variable with control points over a spline - if isinstance(surf_dict["t_over_c"], (int, float)) or surf_dict["t_over_c"].size == 1: - self.set_input_defaults( - f"aero_point.{surf_dict['name']}_perf.t_over_c", val=surf_dict["t_over_c"] * np.ones(ny - 1) - ) - elif surf_dict["t_over_c"].size == ny - 1: - self.set_input_defaults(f"aero_point.{surf_dict['name']}_perf.t_over_c", val=surf_dict["t_over_c"]) - else: - raise ValueError( - f"t_over_c in the surface dict must be either a number or an ndarray " - f"with either one or ny-1 elements, not {surf_dict['t_over_c']}" - ) - - -class PlanformMesh(om.ExplicitComponent): - """ - Generate an OpenAeroStruct mesh based on basic wing design parameters. - Resulting mesh is for a half wing (meant to use with OpenAeroStruct symmetry), - but the input reference area is for the full wing. - - Inputs - ------ - S: float - full planform area (scalar, m^2) - AR: float - aspect ratio (scalar, dimensionless) - taper: float - taper ratio (must be >0 and <=1); tip chord / root chord (scalar, dimensionless) - sweep: float - quarter chord sweep angle (scalar, degrees) - - Outputs - ------- - mesh: ndarray - OpenAeroStruct 3D mesh (num_x x num_y x 3 ndarray, m) - - Options - ------- - num_x: int - number of points in x (streamwise) direction (scalar, dimensionless) - num_y: int - number of points in y (spanwise) direction (scalar, dimensionless) - """ - - def __init__(self, **kwargs): - super().__init__(**kwargs) - self.cite = CITATION - - def initialize(self): - self.options.declare("num_x", default=3, desc="Number of streamwise mesh points") - self.options.declare("num_y", default=7, desc="Number of spanwise (half wing) mesh points") - - def setup(self): - nx = int(self.options["num_x"]) - ny = int(self.options["num_y"]) - - # Generate default mesh - x, y = np.meshgrid(np.linspace(0, 1, nx), np.linspace(-1, 0, ny), indexing="ij") - y *= 5 - mesh = np.zeros((nx, ny, 3)) - mesh[:, :, 0] = x - mesh[:, :, 1] = y - - self.add_input("S", val=10, units="m**2") - self.add_input("AR", val=10) - self.add_input("taper", val=1.0) - self.add_input("sweep", val=10, units="deg") - - self.add_output("mesh", val=mesh, shape=(nx, ny, 3), units="m") - - self.declare_partials("mesh", "*") - - def compute(self, inputs, outputs): - S = inputs["S"] - AR = inputs["AR"] - taper = inputs["taper"] - sweep = inputs["sweep"] - nx = int(self.options["num_x"]) - ny = int(self.options["num_y"]) - - # Compute absolute dimensions from wing geometry spec - half_span = np.sqrt(AR * S) / 2 - c_root = S / (half_span * (1 + taper)) - - # Create baseline square mesh from 0 to 1 in each direction - x_mesh, y_mesh = np.meshgrid(np.linspace(0, 1, nx), np.linspace(-1, 0, ny), indexing="ij") - - # Morph the mesh to fit the desired wing shape - x_mesh *= c_root - y_mesh *= half_span # rectangular wing with correct root chord and wingspan - x_mesh *= np.linspace(taper, 1, ny).reshape(1, ny) # taper wing - x_mesh -= np.linspace(c_root * taper, c_root, ny).reshape(1, ny) / 4 # shift to quarter chord at x=0 - x_mesh += np.linspace(half_span, 0, ny).reshape(1, ny) * np.tan( - np.deg2rad(sweep) - ) # sweep wing at quarter chord - - mesh = np.zeros((nx, ny, 3)) - mesh[:, :, 0] = x_mesh - mesh[:, :, 1] = y_mesh - - outputs["mesh"] = mesh - - def compute_partials(self, inputs, J): - S = inputs["S"] - AR = inputs["AR"] - taper = inputs["taper"] - sweep = inputs["sweep"] - nx = int(self.options["num_x"]) - ny = int(self.options["num_y"]) - - # Compute absolute dimensions from wing geometry spec - half_span = np.sqrt(AR * S) / 2 - c_root = S / (half_span * (1 + taper)) - - # Create baseline square mesh from 0 to 1 in each direction - x_mesh, y_mesh = np.meshgrid(np.linspace(0, 1, nx), np.linspace(-1, 0, ny), indexing="ij") - - # Compute derivatives in a way analogous to forward AD - db_dS = AR / (4 * np.sqrt(AR * S)) - db_dAR = S / (4 * np.sqrt(AR * S)) - dcroot_dS = 1 / (half_span * (1 + taper)) - S / (half_span**2 * (1 + taper)) * db_dS - dcroot_dAR = -S / (half_span**2 * (1 + taper)) * db_dAR - dcroot_dtaper = -S / (half_span * (1 + taper) ** 2) - - dy_dS = y_mesh * db_dS - dy_dAR = y_mesh * db_dAR - - dx_dS = x_mesh * np.linspace(taper, 1, ny).reshape(1, ny) * dcroot_dS - dx_dS -= np.linspace(dcroot_dS * taper, dcroot_dS, ny).reshape(1, ny) / 4 - dx_dS += np.linspace(db_dS, 0, ny).reshape(1, ny) * np.tan(np.deg2rad(sweep)) - - dx_dAR = x_mesh * np.linspace(taper, 1, ny).reshape(1, ny) * dcroot_dAR - dx_dAR -= np.linspace(dcroot_dAR * taper, dcroot_dAR, ny).reshape(1, ny) / 4 - dx_dAR += np.linspace(db_dAR, 0, ny).reshape(1, ny) * np.tan(np.deg2rad(sweep)) - - dx_dtaper = ( - x_mesh * c_root * np.linspace(1, 0, ny).reshape(1, ny) - + x_mesh * np.linspace(taper, 1, ny).reshape(1, ny) * dcroot_dtaper + # Promote the mesh from within OpenAeroStruct + self.promotes( + subsys_name="aero_point", + inputs=[(f"{self.surf_dict['name']}.def_mesh", "ac|geom|wing|OAS_mesh")], ) - dx_dtaper -= ( - np.linspace(c_root, 0, ny).reshape(1, ny) / 4 - + np.linspace(dcroot_dtaper * taper, dcroot_dtaper, ny).reshape(1, ny) / 4 + self.promotes( + subsys_name="aero_point", + inputs=[(f"aero_states.{self.surf_dict['name']}_def_mesh", "ac|geom|wing|OAS_mesh")], ) - dx_dsweep = ( - 0 * x_mesh + np.linspace(half_span, 0, ny).reshape(1, ny) / np.cos(np.deg2rad(sweep)) ** 2 * np.pi / 180.0 - ) - - J["mesh", "S"] = np.dstack((dx_dS, dy_dS, np.zeros((nx, ny)))).flatten() - J["mesh", "AR"] = np.dstack((dx_dAR, dy_dAR, np.zeros((nx, ny)))).flatten() - J["mesh", "taper"] = np.dstack((dx_dtaper, np.zeros((nx, ny)), np.zeros((nx, ny)))).flatten() - J["mesh", "sweep"] = np.dstack((dx_dsweep, np.zeros((nx, ny)), np.zeros((nx, ny)))).flatten() + # Set input defaults for inputs that go to multiple locations + self.set_input_defaults("fltcond|M", 0.1) + self.set_input_defaults("fltcond|alpha", 0.0, units="deg") + self.set_input_defaults("ac|geom|wing|OAS_mesh", dummy_mesh, units="m") + self.set_input_defaults("ac|geom|wing|toverc", np.full(ny - 1, 0.12)) # 12% airfoil thickness by default # Example usage of the drag polar that compares @@ -933,8 +873,8 @@ def compute_partials(self, inputs, J): def example_usage(): # Define parameters nn = 1 - num_x = 3 - num_y = 5 + num_x = 2 + num_y = 4 S = 427.8 # m^2 AR = 9.82 taper = 0.149 @@ -988,17 +928,12 @@ def example_usage(): # Call OpenAeroStruct at the same flight condition to compare prob = om.Problem() - prob.model.add_subsystem("model", VLM(num_x=num_x, num_y=num_y, num_twist=n_twist), promotes=["*"]) + prob.model.add_subsystem("model", VLM(num_x=num_x, num_y=num_y), promotes=["*"]) prob.setup() - # Set values - # Geometry - prob.set_val("ac|geom|wing|S_ref", S, units="m**2") - prob.set_val("ac|geom|wing|AR", AR) - prob.set_val("ac|geom|wing|taper", taper) - prob.set_val("ac|geom|wing|c4sweep", sweep, units="deg") - prob.set_val("ac|geom|wing|twist", twist, units="deg") + # Set mesh value + prob.set_val("ac|geom|wing|OAS_mesh", p.get_val("twisted_mesh", units="m"), units="m") # Flight condition prob.set_val("fltcond|M", M) diff --git a/openconcept/aerodynamics/openaerostruct/mesh_gen.py b/openconcept/aerodynamics/openaerostruct/mesh_gen.py new file mode 100644 index 00000000..6d124a27 --- /dev/null +++ b/openconcept/aerodynamics/openaerostruct/mesh_gen.py @@ -0,0 +1,707 @@ +from copy import deepcopy +import numpy as np +import openmdao.api as om + + +class TrapezoidalPlanformMesh(om.ExplicitComponent): + """ + Generate an OpenAeroStruct mesh based on basic wing design parameters. + Resulting mesh is for a half wing (meant to use with OpenAeroStruct symmetry), + but the input reference area is for the full wing. + + Inputs + ------ + S: float + full planform area (scalar, m^2) + AR: float + aspect ratio (scalar, dimensionless) + taper: float + taper ratio (must be >0 and <=1); tip chord / root chord (scalar, dimensionless) + sweep: float + quarter chord sweep angle (scalar, degrees) + + Outputs + ------- + mesh: ndarray + OpenAeroStruct 3D mesh (num_x + 1 x num_y + 1 x 3 ndarray, m) + + Options + ------- + num_x: int + number of panels in x (streamwise) direction (scalar, dimensionless) + num_y: int + number of panels in y (spanwise) direction (scalar, dimensionless) + """ + + def initialize(self): + self.options.declare("num_x", default=2, desc="Number of streamwise mesh panels") + self.options.declare("num_y", default=6, desc="Number of spanwise (half wing) mesh panels") + + def setup(self): + # Number of coordinates is one more than the number of panels + self.nx = int(self.options["num_x"]) + 1 + self.ny = int(self.options["num_y"]) + 1 + + # Generate default mesh + x, y = np.meshgrid(np.linspace(0, 1, self.nx), np.linspace(-1, 0, self.ny), indexing="ij") + y *= 5 + mesh = np.zeros((self.nx, self.ny, 3)) + mesh[:, :, 0] = x + mesh[:, :, 1] = y + + self.add_input("S", val=10, units="m**2") + self.add_input("AR", val=10) + self.add_input("taper", val=1.0) + self.add_input("sweep", val=10, units="deg") + + self.add_output("mesh", val=mesh, shape=(self.nx, self.ny, 3), units="m") + + self.declare_partials("mesh", "*") + + def compute(self, inputs, outputs): + S = inputs["S"] + AR = inputs["AR"] + taper = inputs["taper"] + sweep = inputs["sweep"] + + # Compute absolute dimensions from wing geometry spec + half_span = np.sqrt(AR * S) / 2 + c_root = S / (half_span * (1 + taper)) + + # Create baseline square mesh from 0 to 1 in each direction + x_mesh, y_mesh = np.meshgrid(np.linspace(0, 1, self.nx), np.linspace(-1, 0, self.ny), indexing="ij") + + # Morph the mesh to fit the desired wing shape + x_mesh *= c_root + y_mesh *= half_span # rectangular wing with correct root chord and wingspan + x_mesh *= np.linspace(taper, 1, self.ny).reshape(1, self.ny) # taper wing + x_mesh -= np.linspace(c_root * taper, c_root, self.ny).reshape(1, self.ny) / 4 # shift to quarter chord at x=0 + x_mesh += np.linspace(half_span, 0, self.ny).reshape(1, self.ny) * np.tan( + np.deg2rad(sweep) + ) # sweep wing at quarter chord + + mesh = np.zeros((self.nx, self.ny, 3)) + mesh[:, :, 0] = x_mesh + mesh[:, :, 1] = y_mesh + + outputs["mesh"] = mesh + + def compute_partials(self, inputs, J): + S = inputs["S"] + AR = inputs["AR"] + taper = inputs["taper"] + sweep = inputs["sweep"] + + # Compute absolute dimensions from wing geometry spec + half_span = np.sqrt(AR * S) / 2 + c_root = S / (half_span * (1 + taper)) + + # Create baseline square mesh from 0 to 1 in each direction + x_mesh, y_mesh = np.meshgrid(np.linspace(0, 1, self.nx), np.linspace(-1, 0, self.ny), indexing="ij") + + # Compute derivatives in a way analogous to forward AD + db_dS = AR / (4 * np.sqrt(AR * S)) + db_dAR = S / (4 * np.sqrt(AR * S)) + dcroot_dS = 1 / (half_span * (1 + taper)) - S / (half_span**2 * (1 + taper)) * db_dS + dcroot_dAR = -S / (half_span**2 * (1 + taper)) * db_dAR + dcroot_dtaper = -S / (half_span * (1 + taper) ** 2) + + dy_dS = y_mesh * db_dS + dy_dAR = y_mesh * db_dAR + + dx_dS = x_mesh * np.linspace(taper, 1, self.ny).reshape(1, self.ny) * dcroot_dS + dx_dS -= np.linspace(dcroot_dS * taper, dcroot_dS, self.ny).reshape(1, self.ny) / 4 + dx_dS += np.linspace(db_dS, 0, self.ny).reshape(1, self.ny) * np.tan(np.deg2rad(sweep)) + + dx_dAR = x_mesh * np.linspace(taper, 1, self.ny).reshape(1, self.ny) * dcroot_dAR + dx_dAR -= np.linspace(dcroot_dAR * taper, dcroot_dAR, self.ny).reshape(1, self.ny) / 4 + dx_dAR += np.linspace(db_dAR, 0, self.ny).reshape(1, self.ny) * np.tan(np.deg2rad(sweep)) + + dx_dtaper = ( + x_mesh * c_root * np.linspace(1, 0, self.ny).reshape(1, self.ny) + + x_mesh * np.linspace(taper, 1, self.ny).reshape(1, self.ny) * dcroot_dtaper + ) + dx_dtaper -= ( + np.linspace(c_root, 0, self.ny).reshape(1, self.ny) / 4 + + np.linspace(dcroot_dtaper * taper, dcroot_dtaper, self.ny).reshape(1, self.ny) / 4 + ) + + dx_dsweep = ( + 0 * x_mesh + + np.linspace(half_span, 0, self.ny).reshape(1, self.ny) / np.cos(np.deg2rad(sweep)) ** 2 * np.pi / 180.0 + ) + + J["mesh", "S"] = np.dstack((dx_dS, dy_dS, np.zeros((self.nx, self.ny)))).flatten() + J["mesh", "AR"] = np.dstack((dx_dAR, dy_dAR, np.zeros((self.nx, self.ny)))).flatten() + J["mesh", "taper"] = np.dstack( + (dx_dtaper, np.zeros((self.nx, self.ny)), np.zeros((self.nx, self.ny))) + ).flatten() + J["mesh", "sweep"] = np.dstack( + (dx_dsweep, np.zeros((self.nx, self.ny)), np.zeros((self.nx, self.ny))) + ).flatten() + + +class SectionPlanformMesh(om.ExplicitComponent): + """ + Generate an OpenAeroStruct mesh based on defined section parameters, similar + to how AVL accepts planforms. The resulting mesh is for a half wing (meant + to use with OpenAeroStruct symmetry), but the input reference area is for the + full wing. + + The sectional properties (x offset, y value, and chord) are nondimensional, + but their relative magnitudes reflect the actual shape. The full planform + is scaled to match the requested planform area. Optionally, this scaling + can be turned off (by setting scale_area=False) and the planform area is + then computed as an output. + + The mesh points are cosine spaced in the chordwise direction and in the + spanwise direction within each trapezoidal region of the wing defined by + pairs of section properties. + + .. code-block:: text + + ----> +y + | + | + v +x + + _ Section 1 Section 2 _ + | | _-' | + | | _-' | + | | _-' | + | _-' chord_sec + x_LE_sec _-' | + | _-' | + | _-' _-------------- _| + | _-' _-' + |_ _-' _-' |--- y_sec ---| + | _-' + |_-' + + Inputs + ------ + S: float + Full planform area; only an input when scale_area is True (the default) (scalar, m^2) + x_LE_sec : float + Streamwise offset of the section's leading edge, starting with the outboard + section (wing tip) and moving inboard toward the root (vector of length + num_sections, m) + y_sec : float + Spanwise location of each section, starting with the outboard section (wing + tip) at the MOST NEGATIVE y value and moving inboard (increasing y value) + toward the root; the user does not provide a value for the root because it + is always 0.0 (vector of length num_sections - 1, m) + chord_sec : float + Chord of each section, starting with the outboard section (wing tip) and + moving inboard toward the root (vector of length num_sections, m) + + Outputs + ------- + mesh: ndarray + OpenAeroStruct 3D mesh (num_x + 1 x sum(num_y) + 1 x 3 ndarray, m) + S: float + Full planform area; only an output when scale_area is False (scalar, m^2) + + Options + ------- + num_x: int + Number of panels in x (streamwise) direction (scalar, dimensionless) + num_y: int or iterable of ints + Number of spanwise panels in the trapezoidal regions between each pair of + adjacent sections; can be specified either as a single integer, where that + same value is used for each region, or as an iterable of integers of length + num_sections - 1 to enable different numbers of spanwise coordinates for each + region (scalar or vector, dimensionless) + num_sections : int + Number of spanwise sections to define planform shape (scalar, dimensionless) + scale_area : bool + Scale the mesh to match a planform area provided as an input + """ # noqa + + def initialize(self): + self.options.declare("num_x", default=3, types=int, desc="Number of streamwise mesh points") + self.options.declare( + "num_y", + default=7, + types=(int, list, tuple, np.ndarray), + desc="Number of spanwise mesh points per trapezoidal region", + ) + self.options.declare( + "num_sections", default=2, types=int, desc="Number of sections along the half span to define" + ) + self.options.declare("scale_area", default=True, types=bool, desc="Scale the planform area") + + def setup(self): + # Process mesh resolution options + self.nx = self.options["num_x"] + 1 # nx is now number of coordinates, not panels + self.n_sec = self.options["num_sections"] + if self.n_sec < 2: + raise ValueError("Must define at least two sections along the half span") + num_y = self.options["num_y"] + if isinstance(num_y, int): + # If it's an int, duplicate the number of sections for each region + self.ny = [num_y] * (self.n_sec - 1) + elif isinstance(num_y, (list, tuple, np.ndarray)): + # If it is an iterable, make sure it's the right length + if len(num_y) != self.n_sec - 1: + raise ValueError("If specified as an iterable, the num_y option must have a length of num_sections") + self.ny = [int(x) for x in num_y] # cast anything that's not an integer to an integer + + self.scale = self.options["scale_area"] + + self.add_input("x_LE_sec", val=0, shape=(self.n_sec,), units="m") + self.add_input("y_sec", val=np.linspace(-5, 0, self.n_sec)[:-1], shape=(self.n_sec - 1,), units="m") + self.add_input("chord_sec", val=1, shape=(self.n_sec,), units="m") + if self.scale: + self.add_input("S", val=10, units="m**2") + self.declare_partials("mesh", "S") + else: + self.add_output("S", val=10, units="m**2") + self.declare_partials("S", ["y_sec", "chord_sec"]) + + # Generate default mesh + self.ny_tot = np.sum(self.ny) + 1 + x, y = np.meshgrid(cos_space(0, 1, self.nx), cos_space(-5, 0, self.ny_tot), indexing="ij") + y *= 5 + mesh = np.zeros((self.nx, self.ny_tot, 3)) + mesh[:, :, 0] = x + mesh[:, :, 1] = y + + self.add_output("mesh", val=mesh, shape=(self.nx, self.ny_tot, 3), units="m") + + # Inputs cache to know if the mesh has been updated or not + self.inputs_cache = None + self.dummy_outputs = {"mesh": np.zeros((self.nx, self.ny_tot, 3))} + + # Compute some values for partial derivatives that don't depend on input values + self.dmesh_dysec_no_S_influence = np.zeros((self.nx * self.ny_tot * 3, self.n_sec - 1), dtype=float) + self.dmesh_dcsec_no_S_influence = np.zeros((self.nx * self.ny_tot * 3, self.n_sec), dtype=float) + self.dmesh_dxLEsec_no_S_influence = np.zeros(self.ny_tot * self.nx * 2, dtype=float) + dmesh_dxLEsec_rows = np.zeros(self.ny_tot * self.nx * 2, dtype=int) # rows for partial derivative sparsity + dmesh_dxLEsec_cols = np.zeros(self.ny_tot * self.nx * 2, dtype=int) # cols for partial derivative sparsity + + # Indices in flattened array + idx_x = 3 * np.arange(self.nx * self.ny_tot) + idx_y = idx_x + 1 + + # Iterate through the defined trapezoidal regions between sections + y_prev = 0 + idx_dmesh_dxLE_prev = 0 + for i_sec in range(self.n_sec - 1): + ny = self.ny[i_sec] + 1 # number of coordinates in the current region (including at the ends) + x_mesh = np.repeat(cos_space(0, 1, self.nx), ny).reshape(self.nx, ny) + + # Derivatives of this trapezoidal region + cos_deriv_start = cos_space_deriv_start(ny) + cos_deriv_end = cos_space_deriv_end(ny) + dymesh_dysec = np.tile(cos_deriv_start, (self.nx, 1)) + dymesh_dysecnext = np.tile(cos_deriv_end, (self.nx, 1)) + dxmesh_dxsec = dymesh_dysec + dxmesh_dxsecnext = dymesh_dysecnext + dxmesh_dcsec = x_mesh * cos_deriv_start + dxmesh_dcsecnext = x_mesh * cos_deriv_end + + # We must be careful not to double count the derivatives at intermediate sections, + # since they are both a beginning and ending point of a trapezoidal region. To do this + # only include the furthest outboard section of the region if this trapezoidal region + # includes the wing tip + y_idx_start = 0 if i_sec == 0 else 1 + # Indices in the mesh (not x, y, z values) corresponding to the coordinates in the current region + idx_mesh = ( + np.tile(np.arange(y_prev + y_idx_start, y_prev + ny), (self.nx, 1)).T + np.arange(self.nx) * self.ny_tot + ).T.flatten() + # No derivative w.r.t. the y value at the root because it's always zero + if i_sec < self.n_sec - 2: + self.dmesh_dysec_no_S_influence[idx_y[idx_mesh], i_sec + 1] += dymesh_dysecnext[ + :, y_idx_start: + ].flatten() + self.dmesh_dysec_no_S_influence[idx_y[idx_mesh], i_sec] += dymesh_dysec[:, y_idx_start:].flatten() + self.dmesh_dcsec_no_S_influence[idx_x[idx_mesh], i_sec + 1] += dxmesh_dcsecnext[:, y_idx_start:].flatten() + self.dmesh_dcsec_no_S_influence[idx_x[idx_mesh], i_sec] += dxmesh_dcsec[:, y_idx_start:].flatten() + + # Derivatives w.r.t. x_LE_sec are sparse, so track rows and cols + idx_end = idx_dmesh_dxLE_prev + idx_mesh.size + dmesh_dxLEsec_rows[idx_dmesh_dxLE_prev:idx_end] = idx_x[idx_mesh] + dmesh_dxLEsec_cols[idx_dmesh_dxLE_prev:idx_end] = i_sec + 1 + self.dmesh_dxLEsec_no_S_influence[idx_dmesh_dxLE_prev:idx_end] = dxmesh_dxsecnext[:, y_idx_start:].flatten() + idx_dmesh_dxLE_prev += idx_mesh.size + + idx_end += idx_mesh.size + dmesh_dxLEsec_rows[idx_dmesh_dxLE_prev:idx_end] = idx_x[idx_mesh] + dmesh_dxLEsec_cols[idx_dmesh_dxLE_prev:idx_end] = i_sec + self.dmesh_dxLEsec_no_S_influence[idx_dmesh_dxLE_prev:idx_end] = dxmesh_dxsec[:, y_idx_start:].flatten() + idx_dmesh_dxLE_prev += idx_mesh.size + + y_prev += ny - 1 + + # Use sparsity for w.r.t. x_LE_sec but don't use sparsity for w.r.t. y_sec and chord_sec since + # they're ~2/3 nonzero regardless of mesh size; still might be better to make y_sec and chord_sec + # zero since derivatives of all z values are zero and the Jacobian will end up as a sparse data + # structure under the hood + self.declare_partials("mesh", ["y_sec", "chord_sec"]) + self.declare_partials("mesh", "x_LE_sec", rows=dmesh_dxLEsec_rows, cols=dmesh_dxLEsec_cols) + + def compute(self, inputs, outputs): + x_sec = inputs["x_LE_sec"] + y_sec = np.hstack((inputs["y_sec"], 0)) + c_sec = inputs["chord_sec"] + + self.inputs_cache = deepcopy(dict(inputs)) + + # Iterate through the defined trapezoidal regions between sections + y_prev = 0 + A = 0.0 # nondimensional area of existing mesh + for i_sec in range(self.n_sec - 1): + ny = self.ny[i_sec] + 1 # number of coordinates in the current region (including at the ends) + x_mesh, y_mesh = np.meshgrid( + cos_space(0, 1, self.nx, dtype=x_sec.dtype), + cos_space(y_sec[i_sec], y_sec[i_sec + 1], ny, dtype=x_sec.dtype), + indexing="ij", + ) + x_mesh *= cos_space(c_sec[i_sec], c_sec[i_sec + 1], ny) + x_mesh += cos_space(x_sec[i_sec], x_sec[i_sec + 1], ny) + + # Because y_prev is incremented by ny - 1, rather than by ny, the furthest inboard + # coordinates from this region will be "overwritten" (with the same values) by + # the next region (not particularly important to know, but might be useful for debugging) + outputs["mesh"][:, y_prev : y_prev + ny, 0] = x_mesh + outputs["mesh"][:, y_prev : y_prev + ny, 1] = y_mesh + + # Add the area of this trapezoidal region to the total nondimensional planform area + A += (y_sec[i_sec + 1] - y_sec[i_sec]) * (c_sec[i_sec] + c_sec[i_sec + 1]) * 0.5 + + y_prev += ny - 1 + + # Zero any z coordinates + outputs["mesh"][:, :, 2] = 0.0 + + # Copy the unscaled mesh for use in the derivative calculation + self.unscaled_flattened_mesh = outputs["mesh"].copy().flatten() + + # Scale the mesh by the reference area + A *= 2 # we're only doing a half wing, so double to get total area + if self.scale: + outputs["mesh"] *= (inputs["S"] / A) ** 0.5 + else: + outputs["S"] = A + + def compute_partials(self, inputs, partials): + y_sec = np.hstack((inputs["y_sec"], 0)) + c_sec = inputs["chord_sec"] + + # Make sure self.usncaled_mesh has been updated with the correct inputs + for input_name in inputs.keys(): + if self.inputs_cache is None or np.any(inputs[input_name] != self.inputs_cache[input_name]): + self.compute(inputs, self.dummy_outputs) + break + + # Load in the initial values of the partials that are independent of the inputs; + # this avoids recomputing them every time compute_partials is called + partials["mesh", "y_sec"][:, :] = self.dmesh_dysec_no_S_influence[:, :] + partials["mesh", "x_LE_sec"][:] = self.dmesh_dxLEsec_no_S_influence[:] + partials["mesh", "chord_sec"][:, :] = self.dmesh_dcsec_no_S_influence[:, :] + + # Iterate through the defined trapezoidal regions between sections + A = 0.0 # nondimensional area of existing mesh + dA_dysec = np.zeros(self.n_sec - 1, dtype=float) + dA_dcsec = np.zeros(self.n_sec, dtype=float) + for i_sec in range(self.n_sec - 1): + # Add the area of this trapezoidal region to the total nondimensional planform area + A += (y_sec[i_sec + 1] - y_sec[i_sec]) * (c_sec[i_sec] + c_sec[i_sec + 1]) * 0.5 + dA_dysec[i_sec] += -0.5 * (c_sec[i_sec] + c_sec[i_sec + 1]) + if i_sec < self.n_sec - 2: + dA_dysec[i_sec + 1] += 0.5 * (c_sec[i_sec] + c_sec[i_sec + 1]) + dA_dcsec[i_sec : i_sec + 2] += 0.5 * (y_sec[i_sec + 1] - y_sec[i_sec]) + + # Scale the mesh by the reference area + A *= 2 # we're only doing a half wing, so double to get total area + dA_dysec *= 2 + dA_dcsec *= 2 + if self.scale: + S = inputs["S"].item() + coeff = (S / A) ** 0.5 + for var in ["y_sec", "chord_sec", "x_LE_sec"]: + partials["mesh", var] *= coeff + partials["mesh", "y_sec"] += np.outer( + -0.5 * self.unscaled_flattened_mesh * S**0.5 * A ** (-1.5), dA_dysec + ) + partials["mesh", "chord_sec"] += np.outer( + -0.5 * self.unscaled_flattened_mesh * S**0.5 * A ** (-1.5), dA_dcsec + ) + partials["mesh", "S"] = 0.5 * S ** (-0.5) / A**0.5 * self.unscaled_flattened_mesh + else: + partials["S", "y_sec"] = dA_dysec + partials["S", "chord_sec"] = dA_dcsec + + +class ThicknessChordRatioInterp(om.ExplicitComponent): + """ + Linearly interpolate thickness to chord ratio defined at each section. + Take the average of the values on either side of a panel to determine + each panel's thickness to chord ratio. + + Inputs + ------ + section_toverc : float + Thickness to chord ratio at each defined section, starting at the + tip and moving to the root (vector of length num_sections, dimensionless) + + Outputs + ------- + panel_toverc : float + Thickness to chord ratio at every streamwise strip of panels in + the mesh (vector of length total panels in y, dimensionless) + + Options + ------- + num_y: int or iterable of ints + Number of spanwise panels in the trapezoidal regions between each pair of + adjacent sections; can be specified either as a single integer, where that + same value is used for each region, or as an iterable of integers of length + num_sections - 1 to enable different numbers of spanwise coordinates for each + region (scalar or vector, dimensionless) + num_sections : int + Number of spanwise sections to define planform shape (scalar, dimensionless) + cos_spacing : bool + mesh is cosine spaced between defined sections, by default True (should be True + if the mesh is generated by SectionPlanformMesh and False if mesh is generated + by TrapezoidalPlanformMesh) + """ + + def initialize(self): + self.options.declare( + "num_y", + default=7, + types=(int, list, tuple, np.ndarray), + desc="Number of spanwise mesh points per trapezoidal region", + ) + self.options.declare( + "num_sections", default=2, types=int, desc="Number of defined sections along the half span" + ) + self.options.declare("cos_spacing", default=True, types=bool, desc="Mesh is cosine spaced within each region") + + def setup(self): + # Process mesh resolution options + self.n_sec = self.options["num_sections"] + if self.n_sec < 2: + raise ValueError("Must define at least two sections along the half span") + num_y = self.options["num_y"] + if isinstance(num_y, int): + # If it's an int, duplicate the number of sections for each region + self.ny = [num_y] * (self.n_sec - 1) + elif isinstance(num_y, (list, tuple, np.ndarray)): + # If it is an iterable, make sure it's the right length + if len(num_y) != self.n_sec - 1: + raise ValueError("If specified as an iterable, the num_y option must have a length of num_sections") + self.ny = [int(x) for x in num_y] # cast anything that's not an integer to an integer + self.ny_tot = np.sum(self.ny) # total number of panels + + self.add_input("section_toverc", val=0.12, shape=(self.n_sec,)) + self.add_output("panel_toverc", val=0.12, shape=(self.ny_tot,)) + + # Compute the partial derivatives (don't depend on inputs since this is linear) + y_prev = 0 + idx_partial_vec = 0 + dpanel_dsection_rows = np.zeros(self.ny_tot * 2, dtype=float) + dpanel_dsection_cols = np.zeros(self.ny_tot * 2, dtype=float) + dpanel_dsection_vals = np.zeros(self.ny_tot * 2, dtype=float) + for i_sec in range(self.n_sec - 1): + ny = self.ny[i_sec] + # First, compute the linearly-interpolated thickness to chord ratio at each y mesh point position + if self.options["cos_spacing"]: + dnodal_dstart = cos_space_deriv_start(ny + 1) + dnodal_dend = cos_space_deriv_end(ny + 1) + else: + dnodal_dstart = np.linspace(1, 0, ny + 1) + dnodal_dend = np.linspace(0, 1, ny + 1) + + # For each panel's t/c, take the average of the nodal t/c's on the panel's two sides + dpanel_dsection_rows[idx_partial_vec : idx_partial_vec + ny] = np.arange(y_prev, y_prev + ny) + dpanel_dsection_cols[idx_partial_vec : idx_partial_vec + ny] = i_sec + dpanel_dsection_vals[idx_partial_vec : idx_partial_vec + ny] = 0.5 * ( + dnodal_dstart[:-1] + dnodal_dstart[1:] + ) + idx_partial_vec += ny + + dpanel_dsection_rows[idx_partial_vec : idx_partial_vec + ny] = np.arange(y_prev, y_prev + ny) + dpanel_dsection_cols[idx_partial_vec : idx_partial_vec + ny] = i_sec + 1 + dpanel_dsection_vals[idx_partial_vec : idx_partial_vec + ny] = 0.5 * (dnodal_dend[:-1] + dnodal_dend[1:]) + idx_partial_vec += ny + + y_prev += ny + + self.declare_partials( + "panel_toverc", + "section_toverc", + rows=dpanel_dsection_rows, + cols=dpanel_dsection_cols, + val=dpanel_dsection_vals, + ) + + def compute(self, inputs, outputs): + tc_section = inputs["section_toverc"] + + # Loop through each region + y_prev = 0 + for i_sec in range(self.n_sec - 1): + ny = self.ny[i_sec] + # First, compute the linearly-interpolated thickness to chord ratio at each y mesh point position + if self.options["cos_spacing"]: + tc_nodal = cos_space(tc_section[i_sec], tc_section[i_sec + 1], ny + 1, dtype=tc_section.dtype) + else: + tc_nodal = np.linspace(tc_section[i_sec], tc_section[i_sec + 1], ny + 1, dtype=tc_section.dtype) + + # For each panel's t/c, take the average of the nodal t/c's on the panel's two sides + outputs["panel_toverc"][y_prev : y_prev + ny] = 0.5 * (tc_nodal[:-1] + tc_nodal[1:]) + + y_prev += ny + + +class SectionLinearInterp(om.ExplicitComponent): + """ + Linearly interpolate a property defined at each section and + evaluate the result at spanwise mesh points. + + Inputs + ------ + property_sec : float + The property to be interpolated at each section definition, + defined from tip to root (vector of length num_sections) + + Outputs + ------- + property_node : float + The linearly interpolated value at each y-coordinate (vector of + length number of different y-coordinates) + + Options + ------- + num_y: int or iterable of ints + Number of spanwise panels in the trapezoidal regions between each pair of + adjacent sections; can be specified either as a single integer, where that + same value is used for each region, or as an iterable of integers of length + num_sections - 1 to enable different numbers of spanwise coordinates for each + region (scalar or vector, dimensionless) + num_sections : int + Number of spanwise sections to define planform shape (scalar, dimensionless) + units : str + Units of the interpolated quantity + cos_spacing : bool + mesh is cosine spaced between defined sections, by default True (should be True + if the mesh is generated by SectionPlanformMesh and False if mesh is generated + by TrapezoidalPlanformMesh) + """ + + def initialize(self): + self.options.declare( + "num_y", + default=7, + types=(int, list, tuple, np.ndarray), + desc="Number of spanwise mesh points per trapezoidal region", + ) + self.options.declare( + "num_sections", default=2, types=int, desc="Number of defined sections along the half span" + ) + self.options.declare("units", default=None, desc="Units of interpolated quantity") + self.options.declare("cos_spacing", default=True, types=bool, desc="Mesh is cosine spaced within each region") + + def setup(self): + # Process mesh resolution options + self.n_sec = self.options["num_sections"] + if self.n_sec < 2: + raise ValueError("Must define at least two sections along the half span") + num_y = self.options["num_y"] + if isinstance(num_y, int): + # If it's an int, duplicate the number of sections for each region + self.ny = [num_y] * (self.n_sec - 1) + elif isinstance(num_y, (list, tuple, np.ndarray)): + # If it is an iterable, make sure it's the right length + if len(num_y) != self.n_sec - 1: + raise ValueError("If specified as an iterable, the num_y option must have a length of num_sections") + self.ny = [int(x) for x in num_y] # cast anything that's not an integer to an integer + self.ny_tot = np.sum(self.ny) + 1 # total number of coordinates + + self.add_input("property_sec", val=0.0, shape=(self.n_sec,), units=self.options["units"]) + self.add_output("property_node", val=0.0, shape=(self.ny_tot,), units=self.options["units"]) + + # Compute the partial derivatives (don't depend on inputs since this is linear) + y_prev = 0 + idx_partial_vec = 0 + num_nonzeros = np.sum(np.asarray(self.ny) - 1) * 2 + self.n_sec + dnodal_dsection_rows = np.zeros(num_nonzeros, dtype=float) + dnodal_dsection_cols = np.zeros(num_nonzeros, dtype=float) + dnodal_dsection_vals = np.zeros(num_nonzeros, dtype=float) + for i_sec in range(self.n_sec - 1): + ny = self.ny[i_sec] + # Compute the linearly-interpolated property at each y mesh point position + if self.options["cos_spacing"]: + dnodal_dstart = cos_space_deriv_start(ny + 1) + dnodal_dend = cos_space_deriv_end(ny + 1) + else: + dnodal_dstart = np.linspace(1, 0, ny + 1) + dnodal_dend = np.linspace(0, 1, ny + 1) + + # Partials from the starting and ending section on nodal properties + idx_start = 0 if i_sec == 0 else 1 + dnodal_dsection_rows[idx_partial_vec : idx_partial_vec + ny - idx_start] = np.arange( + y_prev + idx_start, y_prev + ny + ) + dnodal_dsection_cols[idx_partial_vec : idx_partial_vec + ny - idx_start] = i_sec + dnodal_dsection_vals[idx_partial_vec : idx_partial_vec + ny - idx_start] = dnodal_dstart[idx_start:-1] + idx_partial_vec += ny - idx_start + + dnodal_dsection_rows[idx_partial_vec : idx_partial_vec + ny] = np.arange(y_prev + 1, y_prev + ny + 1) + dnodal_dsection_cols[idx_partial_vec : idx_partial_vec + ny] = i_sec + 1 + dnodal_dsection_vals[idx_partial_vec : idx_partial_vec + ny] = dnodal_dend[1:] + idx_partial_vec += ny + + y_prev += ny + + self.declare_partials( + "property_node", + "property_sec", + rows=dnodal_dsection_rows, + cols=dnodal_dsection_cols, + val=dnodal_dsection_vals, + ) + + def compute(self, inputs, outputs): + prop_sec = inputs["property_sec"] + + # Loop through each region + y_prev = 0 + for i_sec in range(self.n_sec - 1): + ny = self.ny[i_sec] # number of panels in this region + + # Compute the linearly-interpolated property at each y mesh point position + if self.options["cos_spacing"]: + outputs["property_node"][y_prev : y_prev + ny + 1] = cos_space( + prop_sec[i_sec], prop_sec[i_sec + 1], ny + 1, dtype=prop_sec.dtype + ) + else: + outputs["property_node"][y_prev : y_prev + ny + 1] = np.linspace( + prop_sec[i_sec], prop_sec[i_sec + 1], ny + 1, dtype=prop_sec.dtype + ) + + y_prev += ny + + +def cos_space(start, end, num, dtype=float): + """ + Cosine spacing between start and end with num points. + """ + return start + 0.5 * (end - start) * (1 - np.cos(np.linspace(0, np.pi, num, dtype=dtype))) + + +def cos_space_deriv_start(num, dtype=float): + """ + Derivative of cosine spacing output w.r.t. the start value. + """ + return 1 - 0.5 * (1 - np.cos(np.linspace(0, np.pi, num, dtype=dtype))) + + +def cos_space_deriv_end(num, dtype=float): + """ + Derivative of cosine spacing output w.r.t. the end value. + """ + return 0.5 * (1 - np.cos(np.linspace(0, np.pi, num, dtype=dtype))) diff --git a/openconcept/aerodynamics/openaerostruct/tests/test_aerostructural.py b/openconcept/aerodynamics/openaerostruct/tests/test_aerostructural.py index a42cf26c..8bcc116c 100644 --- a/openconcept/aerodynamics/openaerostruct/tests/test_aerostructural.py +++ b/openconcept/aerodynamics/openaerostruct/tests/test_aerostructural.py @@ -50,8 +50,8 @@ def test(self): # Generate mesh to pass to OpenAeroStruct mesh = om.Problem( Aerostruct( - num_x=3, - num_y=7, + num_x=2, + num_y=6, num_twist=twist.size, num_toverc=toverc.size, num_skin=t_skin.size, @@ -75,8 +75,8 @@ def test(self): p = om.Problem( AerostructDragPolar( num_nodes=1, - num_x=3, - num_y=7, + num_x=2, + num_y=6, num_twist=twist.size, num_toverc=toverc.size, num_skin=t_skin.size, @@ -144,8 +144,8 @@ def test_surf_options(self): p = om.Problem( AerostructDragPolar( num_nodes=nn, - num_x=3, - num_y=7, + num_x=2, + num_y=6, num_twist=twist.size, num_toverc=toverc.size, num_skin=t_skin.size, @@ -187,8 +187,8 @@ def test_vectorized(self): p = om.Problem( AerostructDragPolar( num_nodes=nn, - num_x=3, - num_y=7, + num_x=2, + num_y=6, num_twist=twist.size, num_toverc=toverc.size, num_skin=t_skin.size, @@ -242,8 +242,8 @@ def test_defaults(self): p.model.add_subsystem( "comp", OASDataGen( - num_x=3, - num_y=7, + num_x=2, + num_y=6, num_twist=twist.size, num_toverc=toverc.size, num_skin=t_skin.size, @@ -317,7 +317,7 @@ def test_different_surf_options(self): class AerostructTestCase(unittest.TestCase): def get_prob(self, surf_dict={}): p = om.Problem( - Aerostruct(num_x=3, num_y=5, num_twist=2, num_toverc=2, num_skin=2, num_spar=2, surf_options=surf_dict) + Aerostruct(num_x=2, num_y=4, num_twist=2, num_toverc=2, num_skin=2, num_spar=2, surf_options=surf_dict) ) p.setup() p.set_val("fltcond|alpha", 3.0, units="deg") @@ -374,7 +374,7 @@ def test_defaults(self): q = 0.5 * 0.55427276 * 264.20682682**2 nn = 3 p = om.Problem( - AerostructDragPolarExact(num_nodes=nn, num_x=3, num_y=5, num_twist=2, num_toverc=2, num_skin=2, num_spar=2) + AerostructDragPolarExact(num_nodes=nn, num_x=2, num_y=4, num_twist=2, num_toverc=2, num_skin=2, num_spar=2) ) p.model.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, atol=1e-8, rtol=1e-10) p.model.linear_solver = om.DirectSolver() diff --git a/openconcept/aerodynamics/openaerostruct/tests/test_drag_polar.py b/openconcept/aerodynamics/openaerostruct/tests/test_drag_polar.py index c1a83919..57fcb2ac 100644 --- a/openconcept/aerodynamics/openaerostruct/tests/test_drag_polar.py +++ b/openconcept/aerodynamics/openaerostruct/tests/test_drag_polar.py @@ -7,11 +7,11 @@ try: from openaerostruct.geometry.geometry_group import Geometry from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openconcept.aerodynamics.openaerostruct import TrapezoidalPlanformMesh from openconcept.aerodynamics.openaerostruct.drag_polar import ( VLMDataGen, VLM, VLMDragPolar, - PlanformMesh, example_usage, ) @@ -34,24 +34,11 @@ def tearDown(self): def test(self): twist = np.array([-1, -0.5, 2]) - # Generate mesh to pass to OpenAeroStruct - mesh = om.Problem(VLM(num_x=3, num_y=5, num_twist=twist.size)) - mesh.setup() - mesh.set_val("ac|geom|wing|S_ref", 100, units="m**2") - mesh.set_val("ac|geom|wing|AR", 10) - mesh.set_val("ac|geom|wing|taper", 0.1) - mesh.set_val("ac|geom|wing|c4sweep", 20, units="deg") - mesh.set_val("ac|geom|wing|twist", twist, units="deg") - mesh.set_val("fltcond|M", 0.45) - mesh.set_val("fltcond|h", 7.5e3, units="m") - mesh.set_val("fltcond|alpha", 2, units="deg") - mesh.run_model() - p = om.Problem( VLMDragPolar( num_nodes=1, - num_x=3, - num_y=5, + num_x=2, + num_y=4, num_twist=twist.size, Mach_train=np.linspace(0.1, 0.8, 3), alpha_train=np.linspace(-11, 15, 3), @@ -67,10 +54,23 @@ def test(self): p.set_val("ac|geom|wing|taper", 0.1) p.set_val("ac|geom|wing|c4sweep", 20, units="deg") p.set_val("ac|geom|wing|twist", twist, units="deg") + p.set_val("ac|geom|wing|toverc", [0.07, 0.12]) p.set_val("ac|aero|CD_nonwing", 0.01) p.set_val("fltcond|q", 5e3, units="Pa") p.set_val("fltcond|M", 0.45) p.set_val("fltcond|h", 7.5e3, units="m") + p.run_model() + + # Generate mesh to pass to OpenAeroStruct + mesh = om.Problem(VLM(num_x=2, num_y=4)) + mesh.setup() + mesh.set_val("ac|geom|wing|OAS_mesh", p.get_val("twisted_mesh")) + mesh.set_val("ac|geom|wing|toverc", p.get_val("t_over_c_interp.panel_toverc")) + mesh.set_val("fltcond|M", 0.45) + mesh.set_val("fltcond|h", 7.5e3, units="m") + mesh.set_val("fltcond|alpha", 2, units="deg") + mesh.run_model() + p.set_val("fltcond|CL", mesh.get_val("fltcond|CL")) p.run_model() @@ -106,8 +106,8 @@ def test_surf_options(self): p = om.Problem( VLMDragPolar( num_nodes=nn, - num_x=3, - num_y=5, + num_x=2, + num_y=4, num_twist=twist.size, Mach_train=np.linspace(0.1, 0.8, 2), alpha_train=np.linspace(-11, 15, 2), @@ -140,8 +140,8 @@ def test_vectorized(self): p = om.Problem( VLMDragPolar( num_nodes=nn, - num_x=3, - num_y=5, + num_x=2, + num_y=4, num_twist=twist.size, Mach_train=np.linspace(0.1, 0.8, 2), alpha_train=np.linspace(-11, 15, 2), @@ -167,6 +167,103 @@ def test_vectorized(self): # Ensure they're all the same assert_near_equal(p.get_val("drag", units="N"), 37615.14285108 * np.ones(nn), tolerance=1e-10) + def test_section_geometry(self): + nn = 1 + p = om.Problem( + VLMDragPolar( + num_nodes=nn, + num_x=1, + num_y=[2, 3], + num_sections=3, + geometry="section", + Mach_train=np.linspace(0.1, 0.11, 2), + alpha_train=np.linspace(1, 1.1, 2), + alt_train=np.linspace(0, 2, 2), + ) + ) + p.model.nonlinear_solver = om.NewtonSolver(solve_subsystems=True) + p.model.linear_solver = om.DirectSolver() + p.setup() + p.set_val("fltcond|TempIncrement", 0, units="degC") + p.set_val("ac|geom|wing|S_ref", 100, units="m**2") + p.set_val("ac|geom|wing|x_LE_sec", [1, 0.2, 0.0]) + p.set_val("ac|geom|wing|y_sec", [-1, 0.6]) + p.set_val("ac|geom|wing|chord_sec", [0.2, 0.3, 0.3]) + p.set_val("ac|geom|wing|twist", [-1, 0, 1], units="deg") + p.set_val("ac|geom|wing|toverc", [0.1, 0.17, 0.2]) + p.set_val("ac|aero|CD_nonwing", 0.01) + p.set_val("fltcond|q", 5e3, units="Pa") + p.set_val("fltcond|M", 0.105) + p.set_val("fltcond|h", 1, units="m") + p.set_val("fltcond|CL", 0.15579806) + p.run_model() + + vlm = om.Problem(VLM(num_x=1, num_y=5)) + vlm.setup() + vlm.set_val("ac|geom|wing|OAS_mesh", p.get_val("twisted_mesh", units="m"), units="m") + vlm.set_val("fltcond|alpha", 1.05, units="deg") + vlm.set_val("ac|geom|wing|toverc", p.get_val("t_over_c_interp.panel_toverc")) + vlm.set_val("fltcond|M", p.get_val("fltcond|M")) + vlm.set_val("fltcond|h", p.get_val("fltcond|h")) + vlm.run_model() + + # Ensure they're all the same + assert_near_equal(vlm.get_val("fltcond|CL"), 0.15579806, tolerance=1e-3) + assert_near_equal(vlm.get_val("fltcond|CD") + 0.01, p.get_val("aero_surrogate.CD"), tolerance=1e-3) + + def test_mesh_geometry_option(self): + nn = 1 + p = om.Problem() + p.model.add_subsystem( + "mesher", + TrapezoidalPlanformMesh(num_x=1, num_y=2), + promotes=["*"], + ) + p.model.add_subsystem( + "vlm", + VLMDragPolar( + num_nodes=nn, + num_x=1, + num_y=2, + geometry="mesh", + num_twist=3, + Mach_train=np.linspace(0.1, 0.11, 2), + alpha_train=np.linspace(1, 1.1, 2), + alt_train=np.linspace(0, 2, 2), + ), + promotes=["*"], + ) + p.model.connect("mesh", "ac|geom|wing|OAS_mesh") + p.model.nonlinear_solver = om.NewtonSolver(solve_subsystems=True) + p.model.linear_solver = om.DirectSolver() + p.setup() + + p.set_val("fltcond|TempIncrement", 0, units="degC") + p.set_val("S", 100, units="m**2") + p.set_val("AR", 100) + p.set_val("taper", 0.1) + p.set_val("sweep", 20, units="deg") + p.set_val("ac|aero|CD_nonwing", 0.001) + p.set_val("ac|geom|wing|toverc", [0.1, 0.17]) + p.set_val("fltcond|q", 5e3, units="Pa") + p.set_val("fltcond|M", 0.105) + p.set_val("fltcond|h", 1, units="m") + p.set_val("fltcond|CL", 0.10634777) + p.run_model() + + vlm = om.Problem(VLM(num_x=1, num_y=2)) + vlm.setup() + vlm.set_val("ac|geom|wing|OAS_mesh", p.get_val("ac|geom|wing|OAS_mesh", units="m"), units="m") + vlm.set_val("fltcond|alpha", 1.05, units="deg") + vlm.set_val("ac|geom|wing|toverc", p.get_val("ac|geom|wing|toverc")) + vlm.set_val("fltcond|M", p.get_val("fltcond|M")) + vlm.set_val("fltcond|h", p.get_val("fltcond|h")) + vlm.run_model() + + # Ensure they're all the same + assert_near_equal(vlm.get_val("fltcond|CL"), 0.10634777, tolerance=1e-6) + assert_near_equal(vlm.get_val("fltcond|CD") + 0.001, p.get_val("aero_surrogate.CD"), tolerance=1e-4) + @unittest.skipIf(not OAS_installed, "OpenAeroStruct is not installed") class VLMDataGenTestCase(unittest.TestCase): @@ -181,38 +278,41 @@ def tearDown(self): def test_defaults(self): # Regression test - twist = np.array([-1, -0.5, 2]) p = om.Problem() + p.model.add_subsystem( + "mesher", + TrapezoidalPlanformMesh(num_x=1, num_y=2), + promotes=["*"], + ) p.model.add_subsystem( "comp", VLMDataGen( - num_x=3, - num_y=5, - num_twist=twist.size, + num_x=1, + num_y=2, Mach_train=np.linspace(0.1, 0.85, 2), alpha_train=np.linspace(-10, 15, 2), alt_train=np.linspace(0, 15e3, 2), ), promotes=["*"], ) + p.model.connect("mesh", "ac|geom|wing|OAS_mesh") p.setup() p.set_val("fltcond|TempIncrement", 0, units="degC") - p.set_val("ac|geom|wing|S_ref", 100, units="m**2") - p.set_val("ac|geom|wing|AR", 10) - p.set_val("ac|geom|wing|taper", 0.1) - p.set_val("ac|geom|wing|c4sweep", 20, units="deg") - p.set_val("ac|geom|wing|twist", twist, units="deg") + p.set_val("S", 100, units="m**2") + p.set_val("AR", 10) + p.set_val("taper", 0.1) + p.set_val("sweep", 20, units="deg") p.set_val("ac|aero|CD_nonwing", 0.01) p.run_model() CL = np.array( [ - [[-0.79879583, -0.79879583], [1.31170126, 1.31170126]], - [[-0.79879583, -0.79879583], [1.31170126, 1.31170126]], + [[-0.86481567, -0.86481567], [1.28852469, 1.28852469]], + [[-0.86481567, -0.86481567], [1.28852469, 1.28852469]], ] ) CD = np.array( - [[[0.03425776, 0.03647253], [0.06776208, 0.06997685]], [[0.03415836, 0.03597205], [0.20199504, 0.20380873]]] + [[[0.03547695, 0.03770253], [0.05900183, 0.0612274]], [[0.03537478, 0.03719636], [0.18710518, 0.18892676]]] ) assert_near_equal(CL, p.get_val("CL_train"), tolerance=1e-7) @@ -249,25 +349,26 @@ def test_different_surf_options(self): @unittest.skipIf(not OAS_installed, "OpenAeroStruct is not installed") class VLMTestCase(unittest.TestCase): def test_defaults(self): - twist = np.array([-1, -0.5, 2]) - p = om.Problem(VLM(num_x=3, num_y=5, num_twist=twist.size)) + p = om.Problem() + p.model.add_subsystem("mesh", TrapezoidalPlanformMesh(num_x=2, num_y=4), promotes=["*"]) + p.model.add_subsystem("vlm", VLM(num_x=2, num_y=4), promotes=["*"]) + p.model.connect("mesh", "ac|geom|wing|OAS_mesh") p.setup() p.set_val("fltcond|alpha", 2, units="deg") p.set_val("fltcond|M", 0.6) p.set_val("fltcond|h", 5e3, units="m") p.set_val("fltcond|TempIncrement", 0, units="degC") - p.set_val("ac|geom|wing|S_ref", 100, units="m**2") - p.set_val("ac|geom|wing|AR", 10) - p.set_val("ac|geom|wing|taper", 0.1) - p.set_val("ac|geom|wing|c4sweep", 20, units="deg") - p.set_val("ac|geom|wing|twist", twist, units="deg") + p.set_val("S", 100, units="m**2") + p.set_val("AR", 10) + p.set_val("taper", 0.1) + p.set_val("sweep", 20, units="deg") p.run_model() # Run OpenAeroStruct with the same inputs inputs = {} - inputs["mesh"] = p.get_val("mesh.mesh", units="m") - inputs["twist"] = twist + inputs["mesh"] = p.get_val("mesh", units="m") + inputs["twist"] = np.zeros(1) inputs["v"] = p.get_val("airspeed.Utrue", units="m/s") inputs["alpha"] = p.get_val("fltcond|alpha", units="deg") inputs["Mach_number"] = p.get_val("fltcond|M") @@ -280,25 +381,26 @@ def test_defaults(self): assert_near_equal(exact["CD"], p.get_val("fltcond|CD")) def test_wave_drag(self): - twist = np.array([-1, -0.5, 2]) - p = om.Problem(VLM(num_x=3, num_y=5, num_twist=twist.size, surf_options={"with_wave": False})) + p = om.Problem() + p.model.add_subsystem("mesh", TrapezoidalPlanformMesh(num_x=2, num_y=4), promotes=["*"]) + p.model.add_subsystem("vlm", VLM(num_x=2, num_y=4, surf_options={"with_wave": False}), promotes=["*"]) + p.model.connect("mesh", "ac|geom|wing|OAS_mesh") p.setup() p.set_val("fltcond|alpha", 2, units="deg") p.set_val("fltcond|M", 0.85) p.set_val("fltcond|h", 5e3, units="m") p.set_val("fltcond|TempIncrement", 0, units="degC") - p.set_val("ac|geom|wing|S_ref", 100, units="m**2") - p.set_val("ac|geom|wing|AR", 10) - p.set_val("ac|geom|wing|taper", 0.1) - p.set_val("ac|geom|wing|c4sweep", 20, units="deg") - p.set_val("ac|geom|wing|twist", twist, units="deg") + p.set_val("S", 100, units="m**2") + p.set_val("AR", 10) + p.set_val("taper", 0.1) + p.set_val("sweep", 20, units="deg") p.run_model() # Run OpenAeroStruct with the same inputs inputs = {} - inputs["mesh"] = p.get_val("mesh.mesh", units="m") - inputs["twist"] = twist + inputs["mesh"] = p.get_val("mesh", units="m") + inputs["twist"] = np.zeros(1) inputs["v"] = p.get_val("airspeed.Utrue", units="m/s") inputs["alpha"] = p.get_val("fltcond|alpha", units="deg") inputs["Mach_number"] = p.get_val("fltcond|M") @@ -311,25 +413,26 @@ def test_wave_drag(self): assert_near_equal(exact["CD"], p.get_val("fltcond|CD")) def test_viscous_drag(self): - twist = np.array([-1, -0.5, 2]) - p = om.Problem(VLM(num_x=3, num_y=5, num_twist=twist.size, surf_options={"with_viscous": False})) + p = om.Problem() + p.model.add_subsystem("mesh", TrapezoidalPlanformMesh(num_x=2, num_y=4), promotes=["*"]) + p.model.add_subsystem("vlm", VLM(num_x=2, num_y=4, surf_options={"with_viscous": False}), promotes=["*"]) + p.model.connect("mesh", "ac|geom|wing|OAS_mesh") p.setup() p.set_val("fltcond|alpha", 2, units="deg") p.set_val("fltcond|M", 0.85) p.set_val("fltcond|h", 5e3, units="m") p.set_val("fltcond|TempIncrement", 0, units="degC") - p.set_val("ac|geom|wing|S_ref", 100, units="m**2") - p.set_val("ac|geom|wing|AR", 10) - p.set_val("ac|geom|wing|taper", 0.1) - p.set_val("ac|geom|wing|c4sweep", 20, units="deg") - p.set_val("ac|geom|wing|twist", twist, units="deg") + p.set_val("S", 100, units="m**2") + p.set_val("AR", 10) + p.set_val("taper", 0.1) + p.set_val("sweep", 20, units="deg") p.run_model() # Run OpenAeroStruct with the same inputs inputs = {} - inputs["mesh"] = p.get_val("mesh.mesh", units="m") - inputs["twist"] = twist + inputs["mesh"] = p.get_val("mesh", units="m") + inputs["twist"] = np.zeros(1) inputs["v"] = p.get_val("airspeed.Utrue", units="m/s") inputs["alpha"] = p.get_val("fltcond|alpha", units="deg") inputs["Mach_number"] = p.get_val("fltcond|M") @@ -342,25 +445,27 @@ def test_viscous_drag(self): assert_near_equal(exact["CD"], p.get_val("fltcond|CD")) def test_t_over_c(self): - twist = np.array([-1, -0.5, 2]) - p = om.Problem(VLM(num_x=3, num_y=3, num_twist=twist.size, surf_options={"t_over_c": np.array([0.1, 0.2])})) + p = om.Problem() + p.model.add_subsystem("mesh", TrapezoidalPlanformMesh(num_x=2, num_y=2), promotes=["*"]) + p.model.add_subsystem("vlm", VLM(num_x=2, num_y=2), promotes=["*"]) + p.model.connect("mesh", "ac|geom|wing|OAS_mesh") p.setup() p.set_val("fltcond|alpha", 2, units="deg") p.set_val("fltcond|M", 0.85) p.set_val("fltcond|h", 5e3, units="m") p.set_val("fltcond|TempIncrement", 0, units="degC") - p.set_val("ac|geom|wing|S_ref", 100, units="m**2") - p.set_val("ac|geom|wing|AR", 10) - p.set_val("ac|geom|wing|taper", 0.1) - p.set_val("ac|geom|wing|c4sweep", 20, units="deg") - p.set_val("ac|geom|wing|twist", twist, units="deg") + p.set_val("S", 100, units="m**2") + p.set_val("AR", 10) + p.set_val("taper", 0.1) + p.set_val("sweep", 20, units="deg") + p.set_val("ac|geom|wing|toverc", np.array([0.1, 0.2])) p.run_model() # Run OpenAeroStruct with the same inputs inputs = {} - inputs["mesh"] = p.get_val("mesh.mesh", units="m") - inputs["twist"] = twist + inputs["mesh"] = p.get_val("mesh", units="m") + inputs["twist"] = np.zeros(1) inputs["v"] = p.get_val("airspeed.Utrue", units="m/s") inputs["alpha"] = p.get_val("fltcond|alpha", units="deg") inputs["Mach_number"] = p.get_val("fltcond|M") @@ -373,156 +478,6 @@ def test_t_over_c(self): assert_near_equal(exact["CD"], p.get_val("fltcond|CD")) -@unittest.skipIf(not OAS_installed, "OpenAeroStruct is not installed") -class PlanformMeshTestCase(unittest.TestCase): - def test_easy(self): - nx = 3 - ny = 5 - p = om.Problem() - p.model.add_subsystem("comp", PlanformMesh(num_x=nx, num_y=ny), promotes=["*"]) - p.setup() - p.set_val("S", 2, units="m**2") - p.set_val("AR", 2) - p.set_val("taper", 1.0) - p.set_val("sweep", 0.0, units="deg") - p.run_model() - - mesh = np.zeros((nx, ny, 3)) - mesh[:, :, 0], mesh[:, :, 1] = np.meshgrid(np.linspace(-0.25, 0.75, nx), np.linspace(-1, 0, ny), indexing="ij") - - assert_near_equal(p.get_val("mesh", units="m"), mesh) - - partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) - assert_check_partials(partials) - - def test_S_AR(self): - nx = 3 - ny = 5 - p = om.Problem() - p.model.add_subsystem("comp", PlanformMesh(num_x=nx, num_y=ny), promotes=["*"]) - p.setup(force_alloc_complex=True) - p.set_val("S", 48, units="m**2") - p.set_val("AR", 3) - p.set_val("taper", 1.0) - p.set_val("sweep", 0.0, units="deg") - p.run_model() - - mesh = np.zeros((nx, ny, 3)) - mesh[:, :, 0], mesh[:, :, 1] = np.meshgrid(np.linspace(-1, 3, nx), np.linspace(-6, 0, ny), indexing="ij") - - assert_near_equal(p.get_val("mesh", units="m"), mesh) - - partials = p.check_partials(out_stream=None, form="central") - assert_check_partials(partials) - - def test_taper(self): - nx = 2 - ny = 3 - p = om.Problem() - p.model.add_subsystem("comp", PlanformMesh(num_x=nx, num_y=ny), promotes=["*"]) - p.setup() - p.set_val("S", 1.3, units="m**2") - p.set_val("AR", 4 / 1.3) # pick S and AR for half span and root chord of 1 - p.set_val("taper", 0.3) - p.set_val("sweep", 0.0, units="deg") - p.run_model() - - mesh = np.zeros((nx, ny, 3)) - mesh[:, :, 0] = np.array([[-0.075, -0.1625, -0.25], [0.225, 0.4875, 0.75]]) - mesh[:, :, 1] = np.array([[-1, -0.5, 0], [-1, -0.5, 0]]) - - assert_near_equal(p.get_val("mesh", units="m"), mesh) - - partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) - assert_check_partials(partials) - - def test_sweep(self): - nx = 3 - ny = 3 - p = om.Problem() - p.model.add_subsystem("comp", PlanformMesh(num_x=nx, num_y=ny), promotes=["*"]) - p.setup() - p.set_val("S", 2, units="m**2") - p.set_val("AR", 2) - p.set_val("taper", 1.0) - p.set_val("sweep", 45.0, units="deg") - p.run_model() - - mesh = np.zeros((nx, ny, 3)) - mesh[:, :, 0], mesh[:, :, 1] = np.meshgrid(np.linspace(-0.25, 0.75, nx), np.linspace(-1, 0, ny), indexing="ij") - - mesh[:, 0, 0] += 1 - mesh[:, 1, 0] += 0.5 - - assert_near_equal(p.get_val("mesh", units="m"), mesh) - - partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) - assert_check_partials(partials) - - def test_taper_sweep(self): - nx = 2 - ny = 3 - p = om.Problem() - p.model.add_subsystem("comp", PlanformMesh(num_x=nx, num_y=ny), promotes=["*"]) - p.setup() - p.set_val("S", 1.3, units="m**2") - p.set_val("AR", 4 / 1.3) # pick S and AR for half span and root chord of 1 - p.set_val("taper", 0.3) - p.set_val("sweep", 45.0, units="deg") - p.run_model() - - mesh = np.zeros((nx, ny, 3)) - mesh[:, :, 0] = np.array([[-0.075, -0.1625, -0.25], [0.225, 0.4875, 0.75]]) - mesh[:, :, 1] = np.array([[-1, -0.5, 0], [-1, -0.5, 0]]) - mesh[:, 0, 0] += 1 - mesh[:, 1, 0] += 0.5 - - assert_near_equal(p.get_val("mesh", units="m"), mesh) - - partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) - assert_check_partials(partials) - - def test_777ish_regression(self): - nx = 3 - ny = 4 - p = om.Problem() - p.model.add_subsystem("comp", PlanformMesh(num_x=nx, num_y=ny), promotes=["*"]) - p.setup() - p.set_val("S", 427.8, units="m**2") - p.set_val("AR", 9.82) - p.set_val("taper", 0.149) - p.set_val("sweep", 31.6, units="deg") - p.run_model() - - mesh = np.array( - [ - [ - [19.50929722, -32.40754542, 0.0], - [12.04879827, -21.60503028, 0.0], - [4.58829932, -10.80251514, 0.0], - [-2.87219963, 0.0, 0.0], - ], - [ - [20.36521271, -32.40754542, 0.0], - [14.53420835, -21.60503028, 0.0], - [8.70320399, -10.80251514, 0.0], - [2.87219963, 0.0, 0.0], - ], - [ - [21.2211282, -32.40754542, 0.0], - [17.01961843, -21.60503028, 0.0], - [12.81810866, -10.80251514, 0.0], - [8.61659889, 0.0, 0.0], - ], - ] - ) - - assert_near_equal(p.get_val("mesh", units="m"), mesh, tolerance=1e-10) - - partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) - assert_check_partials(partials, atol=2e-5) - - def run_OAS(inputs, with_viscous=True, with_wave=True, t_over_c=None): """ Runs OpenAeroStruct with flight condition and mesh inputs. diff --git a/openconcept/aerodynamics/openaerostruct/tests/test_mesh_gen.py b/openconcept/aerodynamics/openaerostruct/tests/test_mesh_gen.py new file mode 100644 index 00000000..d8592d7e --- /dev/null +++ b/openconcept/aerodynamics/openaerostruct/tests/test_mesh_gen.py @@ -0,0 +1,468 @@ +import unittest +import numpy as np +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials +import openmdao.api as om +from openconcept.aerodynamics.openaerostruct.mesh_gen import ( + TrapezoidalPlanformMesh, + SectionPlanformMesh, + ThicknessChordRatioInterp, + SectionLinearInterp, + cos_space, + cos_space_deriv_start, + cos_space_deriv_end, +) + + +class TrapezoidalPlanformMeshTestCase(unittest.TestCase): + def test_easy(self): + nx = 3 + ny = 5 + p = om.Problem() + p.model.add_subsystem("comp", TrapezoidalPlanformMesh(num_x=nx - 1, num_y=ny - 1), promotes=["*"]) + p.setup() + p.set_val("S", 2, units="m**2") + p.set_val("AR", 2) + p.set_val("taper", 1.0) + p.set_val("sweep", 0.0, units="deg") + p.run_model() + + mesh = np.zeros((nx, ny, 3)) + mesh[:, :, 0], mesh[:, :, 1] = np.meshgrid(np.linspace(-0.25, 0.75, nx), np.linspace(-1, 0, ny), indexing="ij") + + assert_near_equal(p.get_val("mesh", units="m"), mesh) + + partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_S_AR(self): + nx = 3 + ny = 5 + p = om.Problem() + p.model.add_subsystem("comp", TrapezoidalPlanformMesh(num_x=nx - 1, num_y=ny - 1), promotes=["*"]) + p.setup(force_alloc_complex=True) + p.set_val("S", 48, units="m**2") + p.set_val("AR", 3) + p.set_val("taper", 1.0) + p.set_val("sweep", 0.0, units="deg") + p.run_model() + + mesh = np.zeros((nx, ny, 3)) + mesh[:, :, 0], mesh[:, :, 1] = np.meshgrid(np.linspace(-1, 3, nx), np.linspace(-6, 0, ny), indexing="ij") + + assert_near_equal(p.get_val("mesh", units="m"), mesh) + + partials = p.check_partials(out_stream=None, form="central") + assert_check_partials(partials) + + def test_taper(self): + nx = 2 + ny = 3 + p = om.Problem() + p.model.add_subsystem("comp", TrapezoidalPlanformMesh(num_x=nx - 1, num_y=ny - 1), promotes=["*"]) + p.setup() + p.set_val("S", 1.3, units="m**2") + p.set_val("AR", 4 / 1.3) # pick S and AR for half span and root chord of 1 + p.set_val("taper", 0.3) + p.set_val("sweep", 0.0, units="deg") + p.run_model() + + mesh = np.zeros((nx, ny, 3)) + mesh[:, :, 0] = np.array([[-0.075, -0.1625, -0.25], [0.225, 0.4875, 0.75]]) + mesh[:, :, 1] = np.array([[-1, -0.5, 0], [-1, -0.5, 0]]) + + assert_near_equal(p.get_val("mesh", units="m"), mesh) + + partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_sweep(self): + nx = 3 + ny = 3 + p = om.Problem() + p.model.add_subsystem("comp", TrapezoidalPlanformMesh(num_x=nx - 1, num_y=ny - 1), promotes=["*"]) + p.setup() + p.set_val("S", 2, units="m**2") + p.set_val("AR", 2) + p.set_val("taper", 1.0) + p.set_val("sweep", 45.0, units="deg") + p.run_model() + + mesh = np.zeros((nx, ny, 3)) + mesh[:, :, 0], mesh[:, :, 1] = np.meshgrid(np.linspace(-0.25, 0.75, nx), np.linspace(-1, 0, ny), indexing="ij") + + mesh[:, 0, 0] += 1 + mesh[:, 1, 0] += 0.5 + + assert_near_equal(p.get_val("mesh", units="m"), mesh) + + partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_taper_sweep(self): + nx = 2 + ny = 3 + p = om.Problem() + p.model.add_subsystem("comp", TrapezoidalPlanformMesh(num_x=nx - 1, num_y=ny - 1), promotes=["*"]) + p.setup() + p.set_val("S", 1.3, units="m**2") + p.set_val("AR", 4 / 1.3) # pick S and AR for half span and root chord of 1 + p.set_val("taper", 0.3) + p.set_val("sweep", 45.0, units="deg") + p.run_model() + + mesh = np.zeros((nx, ny, 3)) + mesh[:, :, 0] = np.array([[-0.075, -0.1625, -0.25], [0.225, 0.4875, 0.75]]) + mesh[:, :, 1] = np.array([[-1, -0.5, 0], [-1, -0.5, 0]]) + mesh[:, 0, 0] += 1 + mesh[:, 1, 0] += 0.5 + + assert_near_equal(p.get_val("mesh", units="m"), mesh) + + partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_777ish_regression(self): + nx = 3 + ny = 4 + p = om.Problem() + p.model.add_subsystem("comp", TrapezoidalPlanformMesh(num_x=nx - 1, num_y=ny - 1), promotes=["*"]) + p.setup() + p.set_val("S", 427.8, units="m**2") + p.set_val("AR", 9.82) + p.set_val("taper", 0.149) + p.set_val("sweep", 31.6, units="deg") + p.run_model() + + mesh = np.array( + [ + [ + [19.50929722, -32.40754542, 0.0], + [12.04879827, -21.60503028, 0.0], + [4.58829932, -10.80251514, 0.0], + [-2.87219963, 0.0, 0.0], + ], + [ + [20.36521271, -32.40754542, 0.0], + [14.53420835, -21.60503028, 0.0], + [8.70320399, -10.80251514, 0.0], + [2.87219963, 0.0, 0.0], + ], + [ + [21.2211282, -32.40754542, 0.0], + [17.01961843, -21.60503028, 0.0], + [12.81810866, -10.80251514, 0.0], + [8.61659889, 0.0, 0.0], + ], + ] + ) + + assert_near_equal(p.get_val("mesh", units="m"), mesh, tolerance=1e-10) + + partials = p.check_partials(out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials, atol=2e-5) + + +class SectionPlanformMeshTestCase(unittest.TestCase): + def test_hershey_bar(self): + """ + A simple rectangular wing with a span of two and chord of one. + """ + nx = 2 + ny = 2 + x_mesh, y_mesh = np.meshgrid([0, 0.5, 1.0], [-1, -0.5, 0], indexing="ij") + + p = om.Problem() + p.model.add_subsystem("comp", SectionPlanformMesh(num_x=nx, num_y=ny, num_sections=2), promotes=["*"]) + p.setup(force_alloc_complex=True) + + p.set_val("S", 2.0, units="m**2") + p.set_val("x_LE_sec", 0) + p.set_val("y_sec", -1.0) + p.set_val("chord_sec", 1.0) + + p.run_model() + + assert_near_equal(p.get_val("mesh", units="m")[:, :, 0], x_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 1], y_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 2], np.zeros_like(x_mesh)) + + partials = p.check_partials(method="cs", out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_hershey_bar_two_sections(self): + """ + A simple rectangular wing with a span of two and chord of one, but divided into two sections. + """ + nx = 2 + ny = 2 + x_mesh, y_mesh = np.meshgrid(np.linspace(0, 1, nx + 1), np.linspace(-1, 0, ny * 2 + 1), indexing="ij") + + p = om.Problem() + p.model.add_subsystem("comp", SectionPlanformMesh(num_x=nx, num_y=ny, num_sections=3), promotes=["*"]) + p.setup(force_alloc_complex=True) + + p.set_val("S", 2.0, units="m**2") + p.set_val("x_LE_sec", 0) + p.set_val("y_sec", [-1.0, -0.5]) + p.set_val("chord_sec", 1.0) + + p.run_model() + + assert_near_equal(p.get_val("mesh", units="m")[:, :, 0], x_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 1], y_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 2], np.zeros_like(x_mesh)) + + partials = p.check_partials(method="cs", out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_scale_area_off(self): + p = om.Problem() + p.model.add_subsystem( + "comp", SectionPlanformMesh(num_x=3, num_y=3, num_sections=3, scale_area=False), promotes=["*"] + ) + p.setup(force_alloc_complex=True) + + p.set_val("S", 2.0, units="m**2") + p.set_val("x_LE_sec", [2, 0.2, 0.0], units="m") + p.set_val("y_sec", [-1.0, -0.2], units="m") + p.set_val("chord_sec", [0.2, 0.4, 0.4], units="m") + + p.run_model() + + assert_near_equal(p.get_val("S", units="m**2"), 2 * (0.3 * 0.8 + 0.4 * 0.2)) + + partials = p.check_partials(method="cs", out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_chordwise_cos_spacing(self): + nx = 7 + ny = 2 + x_mesh, y_mesh = np.meshgrid(cos_space(0, 1, nx + 1), [-1, -0.5, 0], indexing="ij") + + p = om.Problem() + p.model.add_subsystem("comp", SectionPlanformMesh(num_x=nx, num_y=ny, num_sections=2), promotes=["*"]) + p.setup(force_alloc_complex=True) + + p.set_val("S", 2.0, units="m**2") + p.set_val("x_LE_sec", 0) + p.set_val("y_sec", -1.0) + p.set_val("chord_sec", 1.0) + + p.run_model() + + assert_near_equal(p.get_val("mesh", units="m")[:, :, 0], x_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 1], y_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 2], np.zeros_like(x_mesh)) + + partials = p.check_partials(method="cs", out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_spanwise_cos_spacing(self): + nx = 2 + ny = 7 + x_mesh, y_mesh = np.meshgrid([0.0, 0.5, 1.0], cos_space(-1, 0, ny + 1), indexing="ij") + + p = om.Problem() + p.model.add_subsystem("comp", SectionPlanformMesh(num_x=nx, num_y=ny, num_sections=2), promotes=["*"]) + p.setup(force_alloc_complex=True) + + p.set_val("S", 2.0, units="m**2") + p.set_val("x_LE_sec", 0) + p.set_val("y_sec", -1.0) + p.set_val("chord_sec", 1.0) + + p.run_model() + + assert_near_equal(p.get_val("mesh", units="m")[:, :, 0], x_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 1], y_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 2], np.zeros_like(x_mesh)) + + partials = p.check_partials(method="cs", out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_area_scaling(self): + nx = 2 + ny = 2 + x_mesh, y_mesh = np.meshgrid([0, 1.5, 3.0], [-3, -1.5, 0], indexing="ij") + + p = om.Problem() + p.model.add_subsystem("comp", SectionPlanformMesh(num_x=nx, num_y=ny, num_sections=2), promotes=["*"]) + p.setup(force_alloc_complex=True) + + p.set_val("S", 18.0, units="m**2") + p.set_val("x_LE_sec", 0) + p.set_val("y_sec", -1.0) + p.set_val("chord_sec", 1.0) + + p.run_model() + + assert_near_equal(p.get_val("mesh", units="m")[:, :, 0], x_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 1], y_mesh) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 2], np.zeros_like(x_mesh)) + + partials = p.check_partials(method="cs", out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + def test_BWBish_regression(self): + nx = 3 + num_sec = 4 + x_pos = np.array([480 + 30, 302 + 50, 255, 0]) + y_pos = np.array([-550, -235, -123]) + chord = np.array([30, 100, 577 - 255, 577.0]) + ny = np.array([3, 3, 3]) + + # fmt: off + x_mesh = np.array([[0.88602799, 0.81740425, 0.68015678, 0.61153304, 0.56940328, 0.48514376, 0.443014 , 0.3322605 , 0.1107535 , 0], + [0.89905781, 0.83803481, 0.7159888 , 0.65496579, 0.6369412 , 0.60089202, 0.58286743, 0.49980231, 0.33367206, 0.25060694], + [0.92511746, 0.87929592, 0.78765282, 0.74183128, 0.77201704, 0.83238855, 0.86257431, 0.83488593, 0.77950918, 0.75182081], + [0.93814728, 0.89992647, 0.82348484, 0.78526402, 0.83955495, 0.94813682, 1.00242775, 1.00242775, 1.00242775, 1.00242775]]) + y_mesh = np.array([[-0.95552038, -0.81870724, -0.54508095, -0.4082678 , -0.35962313, -0.26233378, -0.2136891 , -0.16026683, -0.05342228, 0], + [-0.95552038, -0.81870724, -0.54508095, -0.4082678 , -0.35962313, -0.26233378, -0.2136891 , -0.16026683, -0.05342228, 0], + [-0.95552038, -0.81870724, -0.54508095, -0.4082678 , -0.35962313, -0.26233378, -0.2136891 , -0.16026683, -0.05342228, 0], + [-0.95552038, -0.81870724, -0.54508095, -0.4082678 , -0.35962313, -0.26233378, -0.2136891 , -0.16026683, -0.05342228, 0]]) + # fmt: on + + p = om.Problem() + p.model.add_subsystem("mesh", SectionPlanformMesh(num_x=nx, num_y=ny, num_sections=num_sec), promotes=["*"]) + p.setup(force_alloc_complex=True) + p.set_val("S", 0.6, units="m**2") + p.set_val("x_LE_sec", x_pos) + p.set_val("y_sec", y_pos) + p.set_val("chord_sec", chord) + + p.run_model() + + assert_near_equal(p.get_val("mesh", units="m")[:, :, 0], x_mesh, tolerance=1e-6) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 1], y_mesh, tolerance=1e-6) + assert_near_equal(p.get_val("mesh", units="m")[:, :, 2], np.zeros_like(x_mesh)) + + partials = p.check_partials(method="cs", out_stream=None, compact_print=True, show_only_incorrect=True) + assert_check_partials(partials) + + +class ThicknessChordRatioInterpTestCase(unittest.TestCase): + def test_simple(self): + ny = 4 + p = om.Problem() + p.model.add_subsystem( + "comp", ThicknessChordRatioInterp(num_y=ny, num_sections=2, cos_spacing=False), promotes=["*"] + ) + p.setup(force_alloc_complex=True) + + p.set_val("section_toverc", [1, 3]) + p.run_model() + + assert_near_equal(p.get_val("panel_toverc"), 1 + 2 * np.array([0.125, 0.375, 0.625, 0.875])) + + partials = p.check_partials(method="cs") + assert_check_partials(partials) + + def test_cos_spacing(self): + ny = 4 + p = om.Problem() + p.model.add_subsystem( + "comp", ThicknessChordRatioInterp(num_y=ny, num_sections=2, cos_spacing=True), promotes=["*"] + ) + p.setup(force_alloc_complex=True) + + p.set_val("section_toverc", [3, 1]) + p.run_model() + + cos_spacing = cos_space(0, 1, ny + 1) + panel_toverc = 3 - (cos_spacing[:-1] + cos_spacing[1:]) + assert_near_equal(p.get_val("panel_toverc"), panel_toverc) + + partials = p.check_partials(method="cs") + assert_check_partials(partials) + + def test_multiple_sections(self): + ny = np.array([2, 4, 1]) + n_sec = 4 + p = om.Problem() + p.model.add_subsystem( + "comp", ThicknessChordRatioInterp(num_y=ny, num_sections=n_sec, cos_spacing=True), promotes=["*"] + ) + p.setup(force_alloc_complex=True) + + sec_toverc = [0, 2, -2, 1] + p.set_val("section_toverc", sec_toverc) + p.run_model() + + panel_toverc = [] + for i_sec in range(n_sec - 1): + cos_spacing = cos_space(sec_toverc[i_sec], sec_toverc[i_sec + 1], ny[i_sec] + 1) + panel_toverc += list(0.5 * (cos_spacing[:-1] + cos_spacing[1:])) + assert_near_equal(p.get_val("panel_toverc"), panel_toverc) + + partials = p.check_partials(method="cs") + assert_check_partials(partials) + + +class SectionLinearInterpTestCase(unittest.TestCase): + def test_cos_spacing(self): + ny = np.array([3, 2, 1]) + n_sec = 4 + p = om.Problem() + p.model.add_subsystem( + "comp", SectionLinearInterp(num_y=ny, num_sections=n_sec, units="deg", cos_spacing=True), promotes=["*"] + ) + p.setup(force_alloc_complex=True) + + prop = [0, 2, -2, 1] + p.set_val("property_sec", prop) + p.run_model() + + prop_node = np.hstack((cos_space(0, 2, 4), cos_space(2, -2, 3)[1:], [1])) + assert_near_equal(p.get_val("property_node"), prop_node) + + partials = p.check_partials(method="cs") + assert_check_partials(partials) + + def test_linear_spacing(self): + ny = np.array([3, 2, 1]) + n_sec = 4 + p = om.Problem() + p.model.add_subsystem( + "comp", SectionLinearInterp(num_y=ny, num_sections=n_sec, units="deg", cos_spacing=False), promotes=["*"] + ) + p.setup(force_alloc_complex=True) + + prop = [0, 2, -2, 1] + p.set_val("property_sec", prop) + p.run_model() + + prop_node = np.hstack((np.linspace(0, 2, 4), np.linspace(2, -2, 3)[1:], [1])) + assert_near_equal(p.get_val("property_node"), prop_node) + + partials = p.check_partials(method="cs") + assert_check_partials(partials) + + +class CosSpacingTestCase(unittest.TestCase): + def test_simple(self): + num = 7 + correct = (1 - np.cos(np.linspace(0, np.pi, num))) * 0.5 + assert_near_equal(cos_space(0, 1, num), correct) + + def test_shifted(self): + start = 1 + end = -4 + num = 8 + correct = (1 - np.cos(np.linspace(0, np.pi, num))) * 0.5 # 0-1 range + correct = correct * (end - start) + start # stretch and shift + assert_near_equal(cos_space(start, end, num), correct) + + def test_derivs_simple(self): + start = 2 + end = 6 + num = 5 + step = 1e-200 + deriv_start = np.imag(cos_space(start + step * 1j, end, num, dtype=complex)) / step + deriv_end = np.imag(cos_space(start, end + step * 1j, num, dtype=complex)) / step + + assert_near_equal(cos_space_deriv_start(num), deriv_start) + assert_near_equal(cos_space_deriv_end(num), deriv_end) + + +if __name__ == "__main__": + unittest.main() diff --git a/openconcept/examples/B738_VLM_drag.py b/openconcept/examples/B738_VLM_drag.py index d38f3632..df73471d 100644 --- a/openconcept/examples/B738_VLM_drag.py +++ b/openconcept/examples/B738_VLM_drag.py @@ -2,17 +2,19 @@ This work was the basis of the following paper. Please cite it if you use this for your own publication! -@InProceedings{Adler2022a, - author = {Eytan J. Adler and Joaquim R. R. A. Martins}, - title = {Aerostructural wing design optimization considering full mission analysis}, - booktitle = {AIAA SciTech Forum}, - doi = {10.2514/6.2022-0382}, - month = {January}, - year = {2022} +@article{Adler2022d, + author = {Adler, Eytan J. and Martins, Joaquim R. R. A.}, + doi = {10.2514/1.c037096}, + issn = {1533-3868}, + journal = {Journal of Aircraft}, + month = {December}, + publisher = {American Institute of Aeronautics and Astronautics}, + title = {Efficient Aerostructural Wing Optimization Considering Mission Analysis}, + year = {2022} } Eytan Adler (Jan 2022) -""" +""" # noqa import numpy as np @@ -67,7 +69,7 @@ def setup(self): oas_surf_dict["t_over_c"] = acdata["ac"]["geom"]["wing"]["toverc"]["value"] self.add_subsystem( "drag", - VLMDragPolar(num_nodes=nn, num_x=3, num_y=7, num_twist=3, surf_options=oas_surf_dict), + VLMDragPolar(num_nodes=nn, num_x=2, num_y=6, num_twist=3, surf_options=oas_surf_dict), promotes_inputs=[ "fltcond|CL", "fltcond|M", diff --git a/openconcept/examples/B738_aerostructural.py b/openconcept/examples/B738_aerostructural.py index affe1e62..00f2a288 100644 --- a/openconcept/examples/B738_aerostructural.py +++ b/openconcept/examples/B738_aerostructural.py @@ -2,17 +2,19 @@ This work was the basis of the following paper. Please cite it if you use this for your own publication! -@InProceedings{Adler2022a, - author = {Eytan J. Adler and Joaquim R. R. A. Martins}, - title = {Aerostructural wing design optimization considering full mission analysis}, - booktitle = {AIAA SciTech Forum}, - doi = {10.2514/6.2022-0382}, - month = {January}, - year = {2022} +@article{Adler2022d, + author = {Adler, Eytan J. and Martins, Joaquim R. R. A.}, + doi = {10.2514/1.c037096}, + issn = {1533-3868}, + journal = {Journal of Aircraft}, + month = {December}, + publisher = {American Institute of Aeronautics and Astronautics}, + title = {Efficient Aerostructural Wing Optimization Considering Mission Analysis}, + year = {2022} } Eytan Adler (Jan 2022) -""" +""" # noqa import numpy as np @@ -160,8 +162,8 @@ def setup(self): class B738AnalysisGroup(om.Group): def initialize(self): self.options.declare("num_nodes", default=11, desc="Number of analysis points per flight segment") - self.options.declare("num_x", default=3, desc="Aerostructural chordwise nodes") - self.options.declare("num_y", default=7, desc="Aerostructural halfspan nodes") + self.options.declare("num_x", default=2, desc="Aerostructural chordwise nodes") + self.options.declare("num_y", default=6, desc="Aerostructural halfspan nodes") self.options.declare("num_twist", default=3, desc="Number of twist control points") self.options.declare("num_toverc", default=3, desc="Number of t/c control points") self.options.declare("num_skin", default=3, desc="Number of skin control points") @@ -410,8 +412,8 @@ def show_outputs(prob, plots=True): def run_738_analysis(plots=False): num_nodes = 11 global NUM_X, NUM_Y - NUM_X = 3 - NUM_Y = 7 + NUM_X = 2 + NUM_Y = 6 prob = configure_problem(num_nodes) prob.setup(check=False, mode="fwd") set_values(prob, num_nodes) @@ -430,8 +432,8 @@ def run_738_analysis(plots=False): def run_738_optimization(plots=False): num_nodes = 11 global NUM_X, NUM_Y - NUM_X = 3 - NUM_Y = 7 + NUM_X = 2 + NUM_Y = 6 prob = configure_problem(num_nodes) prob.setup(check=True, mode="fwd") set_values(prob, num_nodes) diff --git a/openconcept/thermal/battery_cooling.py b/openconcept/thermal/battery_cooling.py index 34ecd63a..799aafcb 100644 --- a/openconcept/thermal/battery_cooling.py +++ b/openconcept/thermal/battery_cooling.py @@ -169,6 +169,7 @@ class BandolierCoolingSystem(om.ExplicitComponent): t_channel : float Thickness (width) of the cooling channel in the bandolier (scalar, default 1mm) + Outputs ------- dTdt : float diff --git a/openconcept/thermal/heat_exchanger.py b/openconcept/thermal/heat_exchanger.py index cd261d38..407194b4 100644 --- a/openconcept/thermal/heat_exchanger.py +++ b/openconcept/thermal/heat_exchanger.py @@ -489,6 +489,7 @@ class NusseltFromColburnJ(ExplicitComponent): Hydraulic diameter Nusselt number (vector, dimensionless) Nu_dh_hot : float Hydraulic diameter Nusselt number (vector, dimensionless + Options ------- num_nodes : int diff --git a/openconcept/thermal/heat_pipe.py b/openconcept/thermal/heat_pipe.py index 1d3a7dee..c99c8a8f 100644 --- a/openconcept/thermal/heat_pipe.py +++ b/openconcept/thermal/heat_pipe.py @@ -316,6 +316,7 @@ class HeatPipeVaporTempDrop(ExplicitComponent): ------- num_nodes : int Number of analysis points to run, default 1 (scalar, dimensionless) + Other options shouldn't be adjusted since they're for ammonia and there is a hardcoded curve fit also for ammonia in the compute method """ diff --git a/openconcept/utilities/__init__.py b/openconcept/utilities/__init__.py index bfc6bff8..d024c5dd 100644 --- a/openconcept/utilities/__init__.py +++ b/openconcept/utilities/__init__.py @@ -2,7 +2,7 @@ from .dvlabel import DVLabel from .linearinterp import LinearInterpolator from .selector import SelectorComp -from .visualization import plot_trajectory, plot_trajectory_grid +from .visualization import plot_trajectory, plot_trajectory_grid, plot_OAS_mesh, plot_OAS_force_contours # Math utilities from .math import ( diff --git a/openconcept/utilities/math/add_subtract_comp.py b/openconcept/utilities/math/add_subtract_comp.py index 71f916b6..9060eb9f 100644 --- a/openconcept/utilities/math/add_subtract_comp.py +++ b/openconcept/utilities/math/add_subtract_comp.py @@ -65,7 +65,7 @@ def __init__( each input (element-wise addition) val : float or list or tuple or ndarray The initial value of the variable being added in user-defined units. Default is 1.0. - **kwargs : str + kwargs : str Any other arguments to pass to the addition system (same as add_output method for ExplicitComponent) Examples include units (str or None), desc (str) diff --git a/openconcept/utilities/math/combine_split_comp.py b/openconcept/utilities/math/combine_split_comp.py index a90d32db..c9ca948d 100644 --- a/openconcept/utilities/math/combine_split_comp.py +++ b/openconcept/utilities/math/combine_split_comp.py @@ -43,7 +43,7 @@ def __init__(self, output_name=None, input_names=None, vec_sizes=None, length=1, Default is 1 (i.e. a vector of scalars) val : float or list or tuple or ndarray The initial value of the variable being added in user-defined units. Default is 1.0. - **kwargs : str + kwargs : str Any other arguments to pass to the addition system (same as add_output method for ExplicitComponent) Examples include units (str or None), desc (str) @@ -262,7 +262,7 @@ def __init__(self, output_names=None, input_name=None, vec_sizes=None, length=1, Default is 1 (i.e. a vector of scalars) val : float or list or tuple or ndarray The initial value of the variable being added in user-defined units. Default is 1.0. - **kwargs : str + kwargs : str Any other arguments to pass to the addition system (same as add_output method for ExplicitComponent) Examples include units (str or None), desc (str) diff --git a/openconcept/utilities/math/multiply_divide_comp.py b/openconcept/utilities/math/multiply_divide_comp.py index f70f5480..d8fdb4bc 100644 --- a/openconcept/utilities/math/multiply_divide_comp.py +++ b/openconcept/utilities/math/multiply_divide_comp.py @@ -78,7 +78,7 @@ def __init__( input_units : iterable of str Units for each of the input vectors in order. Output units will be dimensionally consistent. - **kwargs : str + kwargs : str Any other arguments to pass to the addition system (same as add_output method for ExplicitComponent) Examples include units (str or None), desc (str) diff --git a/openconcept/utilities/visualization.py b/openconcept/utilities/visualization.py index 233e1c5c..448fc476 100644 --- a/openconcept/utilities/visualization.py +++ b/openconcept/utilities/visualization.py @@ -1,5 +1,6 @@ try: - from matplotlib import pyplot as plt + import matplotlib.pyplot as plt + from matplotlib.collections import LineCollection except ImportError: # don't want a matplotlib dependency on Travis/Appveyor pass @@ -98,3 +99,183 @@ def plot_trajectory_grid( if savefig is not None: fig.tight_layout() plt.savefig(savefig + "_" + str(file_counter) + ".pdf") + + +def plot_OAS_mesh(OAS_mesh, ax=None, set_xlim=True, turn_off_axis=True): + """ + Plots the wing planform mesh from OpenConcept's OpenAeroStruct interface. + + Parameters + ---------- + OAS_mesh : ndarray + The mesh numpy array pulled out of the aerodynamics model (the output + of the mesh component). + ax : matplotlib axis object (optional) + Axis on which to plot the wingbox. If not specified, this function + creates and returns a figure and axis. + set_xlim : bool (optional) + Set the x limits on the axis to just fit the span of the wing. + turn_off_axis : bool (optional) + Turn off the spines, ticks, etc. + + Returns + ------- + fig, ax : matplotlib figure and axis objects + If ax is not specified, returns the created figure and axis. + """ + # Duplicate the half mesh that is used by OAS + mesh = np.hstack((OAS_mesh, OAS_mesh[:, -2::-1, :] * np.array([1, -1, 1]))) + chord_wing_front = mesh[0, :, 0] + span_wing = mesh[0, :, 1] + chord_wing_back = mesh[-1, :, 0] + span = 2 * np.max(span_wing) + + return_ax = False + if ax is None: + # Figure out the size the plot should be + y_range = abs(np.min(chord_wing_front) - np.max(chord_wing_back)) + x_size = 10 + y_size = y_range / span * x_size + + fig, ax = plt.subplots(figsize=(x_size, y_size)) + return_ax = True + + # Plot wing + ax.fill_between( + span_wing, + -chord_wing_front, + -chord_wing_back, + facecolor="#d5e4f5", + zorder=0, + edgecolor="#919191", + clip_on=False, + ) + + # Plot aerodynamic mesh + x = mesh[:, :, 1] + y = -mesh[:, :, 0] + segs1 = np.stack((x, y), axis=2) + segs2 = segs1.transpose(1, 0, 2) + ax.add_collection(LineCollection(segs1, color="#919191", zorder=2, linewidth=0.3)) + ax.add_collection(LineCollection(segs2, color="#919191", zorder=2, linewidth=0.3)) + + # Set final plot details + ax.set_aspect("equal") + if turn_off_axis: + ax.set_axis_off() + if set_xlim: + ax.set_xlim((np.min(span_wing), np.max(span_wing))) + + if return_ax: + return fig, ax + + +def plot_OAS_force_contours( + OAS_mesh, panel_forces, ax=None, set_xlim=True, turn_off_axis=True, wing="both", force_dir=2, **contourf_kwargs +): + """ + Plots contours of the force per area on the surface of the wing. The units are + the force units of panel_forces divided by the square of the lenght units of OAS_mesh. + + Parameters + ---------- + OAS_mesh : ndarray + The mesh numpy array pulled out of the aerodynamics model (the output + of the mesh component). + panel_forces : ndarray + Force from VLM for each panel (the panel_forces output of the VLM component). + ax : matplotlib axis object (optional) + Axis on which to plot the wingbox. If not specified, this function + creates and returns a figure and axis. + set_xlim : bool (optional) + Set the x limits on the axis to just fit the span of the wing, by default True. + turn_off_axis : bool (optional) + Turn off the spines, ticks, etc., by default True. + wing : str (optional) + Which wing to plot, valid options are "left", "right", or "both", by default both. + force_dir : int (optional) + Force direction of which to plot contours. 0 is x force (drag direction), + 1 is y force (spanwise inward), and 2 is z force (lift direction). + contourf_kwargs + Any keyword arguments to pass to matplotlib's contourf function. + + Returns + ------- + c : matplotlib QuadContourSet + Return value from the call to contourf. + fig, ax : matplotlib figure and axis objects + If ax is not specified, returns the created figure and axis. + """ + if wing not in ["left", "right", "both"]: + raise ValueError(f'"{wing}" is not a valid value for wing, must be "left", "right", or "both"') + if force_dir not in [0, 1, 2]: + raise ValueError("force_dir must be either 0, 1, or 2") + + x_mesh = OAS_mesh[:, :, 0] + y_mesh = OAS_mesh[:, :, 1] + chord_wing_front = x_mesh[0, :] + span_wing = y_mesh[0, :] + chord_wing_back = x_mesh[-1, :] + span = np.max(-span_wing) + if wing == "both": + span *= 2 + + return_ax = False + if ax is None: + # Figure out the size the plot should be + y_range = abs(np.min(chord_wing_front) - np.max(chord_wing_back)) + x_size = 10 + y_size = y_range / span * x_size + + fig, ax = plt.subplots(figsize=(x_size, y_size)) + return_ax = True + + # Compute the panel areas + panel_front_widths = y_mesh[:-1, 1:] - y_mesh[:-1, :-1] + panel_back_widths = y_mesh[1:, 1:] - y_mesh[1:, 1:] + panel_avg_widths = 0.5 * (panel_front_widths + panel_back_widths) + panel_left_side_lengths = x_mesh[1:, :-1] - x_mesh[:-1, :-1] + panel_right_side_lengths = x_mesh[1:, 1:] - x_mesh[:-1, 1:] + panel_avg_lengths = 0.5 * (panel_left_side_lengths + panel_right_side_lengths) + panel_areas = panel_avg_widths * panel_avg_lengths + + # Force per area + panel_pressures = panel_forces[:, :, force_dir] / panel_areas + + # Expand the panel forces to be at individual mesh points + le_wt = 0.75 * 0.5 + te_wt = 0.25 * 0.5 + nodal_pressures = np.zeros_like(x_mesh) + nodal_pressures[:-1, :-1] += panel_pressures * le_wt + nodal_pressures[1:, :-1] += panel_pressures * te_wt + nodal_pressures[1:, 1:] += panel_pressures * te_wt + nodal_pressures[:-1, 1:] += panel_pressures * le_wt + + # Adjust the values at the nodes on the edges + nodal_pressures[1:-1, [0, -1]] *= 2 # wing root and tip (but not corners) + nodal_pressures[0, 1:-1] *= 0.5 / le_wt # wing leading edge (not corners) + nodal_pressures[-1, 1:-1] *= 0.5 / te_wt # wing trailing edge (not corners) + nodal_pressures[0, [0, -1]] *= 1 / le_wt # leading edge corners + nodal_pressures[-1, [0, -1]] *= 1 / te_wt # trailing edge corners + + # Plot the contours + if wing == "left": + c = ax.contourf(y_mesh, -x_mesh, nodal_pressures, **contourf_kwargs) + elif wing == "right": + c = ax.contourf(-y_mesh, -x_mesh, nodal_pressures, **contourf_kwargs) + else: + y_mesh_merge = np.hstack((y_mesh, -y_mesh[:, ::-1])) + x_mesh_merge = np.hstack((-x_mesh, -x_mesh[:, ::-1])) + pressures_merge = np.hstack((nodal_pressures, nodal_pressures[:, ::-1])) + c = ax.contourf(y_mesh_merge, x_mesh_merge, pressures_merge, **contourf_kwargs) + + # Set final plot details + ax.set_aspect("equal") + if turn_off_axis: + ax.set_axis_off() + if set_xlim: + ax.set_xlim((np.min(span_wing), np.max(-span_wing))) + + if return_ax: + return c, fig, ax + return c diff --git a/readme.md b/readme.md index 05d0e655..b735db58 100644 --- a/readme.md +++ b/readme.md @@ -63,7 +63,7 @@ OpenConcept is tested regularly on builds with the oldest and latest supported p Please cite this software by reference to the [conference paper](https://www.researchgate.net/publication/326263660_Development_of_a_Conceptual_Design_Model_for_Aircraft_Electric_Propulsion_with_Efficient_Gradients): -Benjamin J. Brelje and Joaquim R.R.A. Martins, "Development of a Conceptual Design Model for Aircraft Electric Propulsion with Efficient Gradients", 2018 AIAA/IEEE Electric Aircraft Technologies Symposium, AIAA Propulsion and Energy Forum, (AIAA 2018-4979) DOI: 10.2514/6.2018-4979 +Benjamin J. Brelje and Joaquim R. R. A. Martins, "Development of a Conceptual Design Model for Aircraft Electric Propulsion with Efficient Gradients", 2018 AIAA/IEEE Electric Aircraft Technologies Symposium, AIAA Propulsion and Energy Forum, (AIAA 2018-4979) DOI: 10.2514/6.2018-4979 ``` @inproceedings{Brelje2018a, @@ -79,15 +79,17 @@ Benjamin J. Brelje and Joaquim R.R.A. Martins, "Development of a Conceptual Desi If using the integrated OpenAeroStruct VLM or aerostructural aerodynamic models, please cite the following [conference paper](https://www.researchgate.net/publication/357559489_Aerostructural_wing_design_optimization_considering_full_mission_analysis): -Eytan J. Adler and Joaquim R.R.A. Martins, "Aerostructural wing design optimization considering full mission analysis", 2022 AIAA SciTech Forum, San Diego, CA, January 2022. DOI: 10.2514/6.2022-0382 +Eytan J. Adler and Joaquim R. R. A. Martins, "Efficient Aerostructural Wing Optimization Considering Mission Analysis", Journal of Aircraft, 2022. DOI: 10.2514/1.c037096 ``` -@inproceedings{Adler2022a, - author = {Eytan J. Adler and Joaquim R. R. A. Martins}, - title = {Aerostructural wing design optimization considering full mission analysis}, - booktitle = {AIAA SciTech Forum}, - doi = {10.2514/6.2022-0382}, - month = {January}, - year = {2022} +@article{Adler2022d, + author = {Adler, Eytan J. and Martins, Joaquim R. R. A.}, + doi = {10.2514/1.c037096}, + issn = {1533-3868}, + journal = {Journal of Aircraft}, + month = {December}, + publisher = {American Institute of Aeronautics and Astronautics}, + title = {Efficient Aerostructural Wing Optimization Considering Mission Analysis}, + year = {2022} } ```