diff --git a/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb b/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb index 77621f9e2..b44ec609f 100644 --- a/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb +++ b/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb @@ -582,7 +582,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.3" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/aviary/subsystems/mass/simple_mass/__init__.py b/aviary/subsystems/mass/simple_mass/__init__.py new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/__init__.py @@ -0,0 +1 @@ + diff --git a/aviary/subsystems/mass/simple_mass/fuselage.py b/aviary/subsystems/mass/simple_mass/fuselage.py new file mode 100644 index 000000000..3abb31ca0 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/fuselage.py @@ -0,0 +1,184 @@ +from pathlib import Path + +import jax +import jax.numpy as jnp +import jax.scipy.interpolate as jinterp +import numpy as np +import openmdao.api as om +import openmdao.jax as omj + +from aviary.subsystems.mass.simple_mass.materials_database import materials +from aviary.utils.named_values import get_keys +from aviary.variable_info.functions import add_aviary_input, add_aviary_output +from aviary.variable_info.variables import Aircraft + + +class FuselageMass(om.JaxExplicitComponent): + def initialize(self): + self.options.declare('num_sections', types=int, default=10) + + self.options.declare('material', default='Aluminum Oxide', values=list(get_keys(materials))) + + self.options.declare( + 'fuselage_data_file', + types=(Path, str), + default=None, + allow_none=True, + desc='optional data file of fuselage geometry', + ) + + self.options.declare( + 'hollow_fuselage', + default=True, + desc='flag for if the fuselage should be hollow (wall thickness is not used if False)', + ) + + # TODO FunctionType is not defined? + # self.options.declare( + # 'custom_fuselage_function', + # types=FunctionType, + # default=None, + # allow_none=True, + # desc='optional custom function generation for fuselage geometry', + # ) + + def setup(self): + # Inputs + add_aviary_input(self, Aircraft.Fuselage.LENGTH, units='m') + + self.add_input('base_diameter', val=0.4, units='m') # no aviary input + self.add_input('tip_diameter', val=0.2, units='m') # no aviary input + self.add_input('curvature', val=0.0, units='m') # 0 for straight, positive for upward curve + self.add_input('thickness', val=0.05, units='m') # Wall thickness of the fuselage + # Allow for asymmetry in the y and z axes -- this value acts as a slope for linear variation along these axes + self.add_input('y_offset', val=0.0, units='m') + self.add_input('z_offset', val=0.0, units='m') + + # Outputs + add_aviary_output(self, Aircraft.Fuselage.MASS, units='kg') + + def compute_primal( + self, + aircraft__fuselage__length, + base_diameter, + tip_diameter, + curvature, + thickness, + y_offset, + z_offset, + ): + is_hollow = self.options['hollow_fuselage'] + + self.validate_inputs( + aircraft__fuselage__length, base_diameter, tip_diameter, thickness, is_hollow + ) + + custom_fuselage_data_file = self.options['fuselage_data_file'] + material = self.options['material'] + num_sections = self.options['num_sections'] + + density = materials.get_val(material, 'kg/m**3') + + section_locations = jnp.linspace(0, aircraft__fuselage__length, num_sections) + + aircraft__fuselage__mass = 0 + total_moment_x = 0 + total_moment_y = 0 + total_moment_z = 0 + + # Load fuselage data file if present + if custom_fuselage_data_file: + try: + # Load the file + custom_data = np.loadtxt(custom_fuselage_data_file) + except FileNotFoundError as e: + raise FileNotFoundError(f'Fuselage data file {e}') + else: + fuselage_locations = custom_data[:, 0] + fuselage_diameters = custom_data[:, 1] + # TODO: OM interp is much more performant than scipy, use metamodel here + interpolate_diameter = jinterp.RegularGridInterpolator( + fuselage_locations, fuselage_diameters, method='linear' + ) + else: + interpolate_diameter = None + + # Loop through each section + for location in section_locations: + # FunctionType is not defined, so "custom_fuselage_function" is currently broken + # if self.options['custom_fuselage_function'] is not None: + # section_diameter = self.options['custom_fuselage_function'](location) + # should be elif below once fixed + if self.options['fuselage_data_file'] is not None and interpolate_diameter is not None: + section_diameter = interpolate_diameter(location) + else: + section_diameter = ( + base_diameter + + ((tip_diameter - base_diameter) / aircraft__fuselage__length) * location + ) + + outer_radius = section_diameter / 2.0 + inner_radius = jnp.where(is_hollow, omj.smooth_max(0, outer_radius - thickness), 0) + + section_volume = ( + jnp.pi + * (outer_radius**2 - inner_radius**2) + * (aircraft__fuselage__length / num_sections) + ) + section_weight = density * section_volume + + centroid_x = jnp.where(tip_diameter / base_diameter != 1, (3 / 4) * location, location) + centroid_y = y_offset * (1 - location / aircraft__fuselage__length) + centroid_z = ( + z_offset * (1 - location / aircraft__fuselage__length) + + curvature * location**2 / aircraft__fuselage__length + ) + + aircraft__fuselage__mass += section_weight + total_moment_x += centroid_x * section_weight + total_moment_y += centroid_y * section_weight + total_moment_z += centroid_z * section_weight + + return aircraft__fuselage__mass + + def validate_inputs(self, length, base_diameter, tip_diameter, thickness, is_hollow): + jax.lax.cond( + length[0] <= 0, + lambda: jax.debug.callback( + raise_error, om.AnalysisError('Length must be a positive value.') + ), + lambda: None, + ) + jax.lax.cond( + base_diameter[0] <= 0, + lambda: jax.debug.callback( + raise_error, om.AnalysisError('Base Diameter must be a positive value.') + ), + lambda: None, + ) + jax.lax.cond( + tip_diameter[0] <= 0, + lambda: jax.debug.callback( + raise_error, om.AnalysisError('Tip Diameter must be a positive value.') + ), + lambda: None, + ) + jax.lax.cond( + thickness[0] <= 0, + lambda: jax.debug.callback( + raise_error, om.AnalysisError('Thickness must be a positive value.') + ), + lambda: None, + ) + + jax.lax.cond( + is_hollow and thickness[0] >= base_diameter[0] / 2, + lambda: jax.debug.callback( + raise_error, om.AnalysisError('Wall thickness is too large for a hollow fuselage.') + ), + lambda: None, + ) + + +def raise_error(exception): + raise exception diff --git a/aviary/subsystems/mass/simple_mass/mass_builder.py b/aviary/subsystems/mass/simple_mass/mass_builder.py new file mode 100644 index 000000000..c51209823 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/mass_builder.py @@ -0,0 +1,25 @@ +from aviary.subsystems.mass.mass_builder import MassBuilderBase +from aviary.subsystems.mass.simple_mass.mass_premission import SimpleMassPremission + +""" +Define subsystem builder for Aviary simple mass. + +Classes +-------- +SimpleMassBuilderBase: the interface for the simple mass subsystem builder. +""" + +_default_name = 'simple_mass' + + +class SimpleMassBuilder(MassBuilderBase): + """Base mass builder.""" + + def __init__(self, name=None): + if name is None: + name = _default_name + + super().__init__(name=name) + + def build_pre_mission(self, aviary_inputs): + return SimpleMassPremission() diff --git a/aviary/subsystems/mass/simple_mass/mass_premission.py b/aviary/subsystems/mass/simple_mass/mass_premission.py new file mode 100644 index 000000000..fe4285d31 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/mass_premission.py @@ -0,0 +1,30 @@ +import openmdao.api as om + +from aviary.subsystems.mass.simple_mass.fuselage import FuselageMass +from aviary.subsystems.mass.simple_mass.mass_summation import SimpleMassSummation +from aviary.subsystems.mass.simple_mass.tail import HorizontalTailMass, VerticalTailMass +from aviary.subsystems.mass.simple_mass.wing import WingMass + + +class SimpleMassPremission(om.Group): + """ + Pre-mission group of top-level mass estimation groups and components for + the simple small-scale aircraft mass build-up. + """ + + def setup(self): + self.add_subsystem('Wing', WingMass(), promotes_inputs=['*'], promotes_outputs=['*']) + + self.add_subsystem( + 'Fuselage', FuselageMass(), promotes_inputs=['*'], promotes_outputs=['*'] + ) + + self.add_subsystem( + 'HorizontalTail', HorizontalTailMass(), promotes_inputs=['*'], promotes_outputs=['*'] + ) + + self.add_subsystem( + 'VerticalTail', VerticalTailMass(), promotes_inputs=['*'], promotes_outputs=['*'] + ) + + self.add_subsystem('mass_summation', SimpleMassSummation(), promotes=['*']) diff --git a/aviary/subsystems/mass/simple_mass/mass_summation.py b/aviary/subsystems/mass/simple_mass/mass_summation.py new file mode 100644 index 000000000..ada5f98de --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/mass_summation.py @@ -0,0 +1,44 @@ +import openmdao.api as om + +from aviary.variable_info.functions import add_aviary_input +from aviary.variable_info.variables import Aircraft + + +class SimpleMassSummation(om.Group): + """ + + Group to compute various design masses for this mass group. + + This group will be expanded greatly as more subsystems are created. + + """ + + def setup(self): + self.add_subsystem( + 'structure_mass', StructureMass(), promotes_inputs=['*'], promotes_outputs=['*'] + ) + + +class StructureMass(om.ExplicitComponent): + def setup(self): + add_aviary_input(self, Aircraft.Wing.MASS, val=0.0, units='kg') + add_aviary_input(self, Aircraft.Fuselage.MASS, val=0.0, units='kg') + add_aviary_input(self, Aircraft.HorizontalTail.MASS, val=0.0, units='kg') + add_aviary_input(self, Aircraft.VerticalTail.MASS, val=0.0, units='kg') + + # More masses can be added, i.e., tail, spars, flaps, etc. as needed + + self.add_output(Aircraft.Design.STRUCTURE_MASS, val=0.0, units='kg') + + def setup_partials(self): + self.declare_partials(Aircraft.Design.STRUCTURE_MASS, '*', val=1) + + def compute(self, inputs, outputs): + wing_mass = inputs[Aircraft.Wing.MASS] + fuselage_mass = inputs[Aircraft.Fuselage.MASS] + htail_mass = inputs[Aircraft.HorizontalTail.MASS] + vtail_mass = inputs[Aircraft.VerticalTail.MASS] + + structure_mass = wing_mass + fuselage_mass + htail_mass + vtail_mass + + outputs[Aircraft.Design.STRUCTURE_MASS] = structure_mass diff --git a/aviary/subsystems/mass/simple_mass/materials_database.py b/aviary/subsystems/mass/simple_mass/materials_database.py new file mode 100644 index 000000000..c63247f29 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/materials_database.py @@ -0,0 +1,72 @@ +""" +Database for various material densities that are to be used for mass calculations for small aircraft in particular. + +This database will be expanded as needed. + +""" + +from aviary.utils.named_values import NamedValues + +materials = NamedValues() + +""" +All densities below came from https://tpsx.arc.nasa.gov/MaterialsDatabase + +""" + +# Wood +materials.set_val('Balsa', 130, units='kg/m**3') +materials.set_val('Cypress', 460, units='kg/m**3') +materials.set_val('Mahogany', 540, units='kg/m**3') +materials.set_val('Maple', 710, units='kg/m**3') +materials.set_val('Teak', 640, units='kg/m**3') + +# Aluminum Compounds and Alloys +materials.set_val('Aluminum Oxide', 3400, units='kg/m**3') +materials.set_val('2024-T8XX', 2800, units='kg/m**3') # aircraft-grade strength Aluminum alloy +materials.set_val('2219-T8XX', 2810, units='kg/m**3') # Exceptionally strong Aluminum alloy +materials.set_val('2024-T6', 2770, units='kg/m**3') # Another Aluminum alloy +materials.set_val('Aluminum Foam', 1300, units='kg/m**3') + +# Steel +materials.set_val('Stainless Steel 17-4 PH', 7830, units='kg/m**3') # 17-4 PH stainless steel +materials.set_val('Stainless Steel-AISI 302', 8060, units='kg/m**3') # AISI 302 +materials.set_val('Stainless Steel-AISI 304', 7900, units='kg/m**3') # AISI 304 +materials.set_val('Steel Alloy Cast', 7830, units='kg/m**3') # General steel alloy cast +materials.set_val('Steel 321', 8030, units='kg/m**3') # Steel type 321 + +# Carbon Fibers / Carbon - Silicon Fibers +materials.set_val('Carbon/Silicon-Carbide', 2080, units='kg/m**3') # Carbon fiber reinforced SiC +materials.set_val( + 'Silicon-Carbide/Silicon-Carbide', 2400, units='kg/m**3' +) # SiC fiber reinforced SiC matrix +materials.set_val('Advanced Carbon-Carbon Composite', 1600, units='kg/m**3') # ACC +materials.set_val('Reinforced Carbon-Carbon', 1580, units='kg/m**3') +materials.set_val( + 'Reinforced Carbon-Carbon Composite', 1580, units='kg/m**3' +) # Generally, ACC is better, but RCC is slightly cheaper + +""" +Below are miscellaneous values that could be of importance, particularly for small aircraft. + +These values were found from a variety of sources, and depending on the source/brand, the density +could be slightly different. For some cases, temperature of the material also matters (typically +the values are provided as a relative density). If there is a temperature dependence from the source, +it will be noted as a comment next to the line where the material value is set. Below are some sources +for various values. + +The values below were not explicity listed from the above source. + +Wood glue: https://www.gorillatough.com/wp-content/uploads/Gorilla-Wood-Glue-v1.2.pdf + +EPS Foam: https://www.abtfoam.com/wp-content/uploads/2020/05/EPS-Standard-Sheet-Sizes-Densities-and-R-values.pdf + Note that there is a density range given, along with different types. The density value used is for Type I, + and the value given is the average of the minimum and maximum within the range provided. The base unit in + this document is pcf for the density. It was converted to kg/m^3 for the actual value input. + +""" + +materials.set_val( + 'Wood Glue', 1080, units='kg/m**3' +) # Relative density value -- corresponds to 25 C (77 F) +materials.set_val('EPS Foam', 16.3388, units='kg/m**3') diff --git a/aviary/subsystems/mass/simple_mass/tail.py b/aviary/subsystems/mass/simple_mass/tail.py new file mode 100644 index 000000000..70f6b07b6 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/tail.py @@ -0,0 +1,335 @@ +from pathlib import Path + +import jax.numpy as jnp +import numpy as np +import openmdao.api as om + +from aviary.subsystems.mass.simple_mass.materials_database import materials +from aviary.subsystems.mass.simple_mass.utils import ( + airfoil_camber_line, + airfoil_thickness, + extract_airfoil_features, + precompute_airfoil_geometry, +) +from aviary.utils.functions import get_path +from aviary.utils.named_values import get_keys +from aviary.variable_info.functions import add_aviary_input, add_aviary_output +from aviary.variable_info.variables import Aircraft + +try: + from quadax import quadgk +except ImportError: + raise ImportError( + "quadax package not found. You can install it by running 'pip install quadax'." + ) + + +# BUG this component is not currently working - mass is not being properly computed +class HorizontalTailMass(om.JaxExplicitComponent): + def initialize(self): + self.options.declare( + 'NACA_digits', + default='2412', + types=(str, int), + desc='4 digit code for NACA airfoil of tail', + ) + + self.options.declare( + 'airfoil_data_file', + default=None, + types=(str, Path), + desc='File path for airfoil coordinates (overrides NACA digits)', + ) + + self.options.declare( + 'material', default='Balsa', values=list(get_keys(materials)), desc='Material type' + ) + + self.options.declare('num_sections', default=10, desc='Number of sections for enumeration') + + self.camber = 0 + self.camber_location = 0 + self.max_thickness = 0 + self.camber_line = 0 + self.thickness = 0 + + def setup(self): + # Inputs + # TODO most of these are unused?!? + add_aviary_input(self, Aircraft.HorizontalTail.SPAN, units='m', desc='Tail span') + add_aviary_input( + self, Aircraft.HorizontalTail.ROOT_CHORD, units='m', desc='Root chord length' + ) + + self.add_input('tip_chord_tail', val=0.8, units='m', desc='Tip chord length') + self.add_input( + 'thickness_ratio', val=0.12, desc='Max thickness to chord ratio for NACA airfoil' + ) + self.add_input('skin_thickness', val=0.002, units='m', desc='Skin panel thickness') + self.add_input( + 'twist_tail', + val=jnp.zeros(self.options['num_sections']), + units='deg', + desc='Twist distribution', + ) + + # Outputs + add_aviary_output( + self, Aircraft.HorizontalTail.MASS, units='kg', desc='Total mass of the tail' + ) + + # File check + airfoil_file = self.options['airfoil_data_file'] + if airfoil_file is not None: + airfoil_file = get_path(airfoil_file) + airfoil_data = np.loadtxt(airfoil_file, skiprows=1) # Assume a header + + x_coords, y_coords = airfoil_data[:, 0], airfoil_data[:, 1] + + ( + self.camber, + self.camber_location, + self.max_thickness, + self.camber_line, + self.thickness, + ) = extract_airfoil_features(x_coords, y_coords) + + else: + NACA_digits = str(self.options['NACA_digits']) + # Parse the NACA airfoil type (4-digit) + self.camber = int(NACA_digits[0]) / 100.0 # Maximum camber + self.camber_location = int(NACA_digits[1]) / 10.0 # Location of max camber + self.max_thickness = int(NACA_digits[2:4]) / 100.0 # Max thickness + + def get_self_statics(self): + return ( + self.camber, + self.camber_location, + self.max_thickness, + self.camber_line, + self.thickness, + self.options['material'], + self.options['num_sections'], + self.options['NACA_digits'], + ) + + def compute_primal( + self, + aircraft__horizontal_tail__span, + aircraft__horizontal_tail__root_chord, + tip_chord_tail, + thickness_ratio, + skin_thickness, + twist_tail, + ): + material = self.options['material'] + density = materials.get_val(material, 'kg/m**3') + airfoil_file = self.options['airfoil_data_file'] + num_sections = self.options['num_sections'] + camber = self.camber + camber_location = self.camber_location + max_thickness = self.max_thickness + camber_line = self.camber_line + thickness = self.thickness + + # This is just so that the differentiation and unittest do not break. + aircraft__horizontal_tail__mass = 0.0 * thickness_ratio + + # TODO unused?? + span_locations = jnp.linspace(0, aircraft__horizontal_tail__span, num_sections) + + # Get x_points and dx for later + x_points, dx = precompute_airfoil_geometry(num_sections) + + # Thickness distribution + thickness_dist = airfoil_thickness(x_points, max_thickness) + + if airfoil_file is not None: + aircraft__horizontal_tail__mass, _ = quadgk( + density * 2 * thickness * jnp.sqrt(1 + jnp.gradient(camber_line) ** 2), + [0, 1], + epsabs=1e-9, + epsrel=1e-9, + ) + else: + total_mass_first_part, _ = quadgk( + lambda x: density + * 2 + * jnp.atleast_1d(airfoil_thickness(x, max_thickness)) + * jnp.sqrt( + 1 + ((camber / camber_location**2) * (2 * camber_location - 2 * x)) ** 2 + ), + [0, camber_location], + epsabs=1e-9, + epsrel=1e-9, + ) + total_mass_second_part, _ = quadgk( + lambda x: density + * 2 + * jnp.atleast_1d(airfoil_thickness(x, max_thickness)) + * jnp.sqrt( + 1 + (camber / (1 - camber_location) ** 2 * (2 * camber_location - 2 * x)) ** 2 + ), + [camber_location, 1], + epsabs=1e-9, + epsrel=1e-9, + ) + + aircraft__horizontal_tail__mass = total_mass_first_part + total_mass_second_part + + return aircraft__horizontal_tail__mass + + +class VerticalTailMass(om.JaxExplicitComponent): + def initialize(self): + self.options.declare( + 'NACA_digits', + default='2412', + types=(str, int), + desc='4 digit code for NACA airfoil of tail', + ) + + self.options.declare( + 'airfoil_data_file', + default=None, + types=(str, Path), + desc='File path for airfoil coordinates (overrides NACA digits)', + ) + + self.options.declare( + 'material', default='Balsa', values=list(get_keys(materials)), desc='Material type' + ) + + self.options.declare('num_sections', default=10, desc='Number of sections for enumeration') + + self.camber = 0 + self.camber_location = 0 + self.max_thickness = 0 + self.camber_line = 0 + self.thickness = 0 + + def setup(self): + # Inputs + # TODO most of these are unused?!? + add_aviary_input(self, Aircraft.VerticalTail.SPAN, units='m', desc='Tail span') + add_aviary_input( + self, Aircraft.VerticalTail.ROOT_CHORD, units='m', desc='Root chord length' + ) + + self.add_input('tip_chord_tail', val=0.8, units='m', desc='Tip chord length') + self.add_input( + 'thickness_ratio', val=0.12, desc='Max thickness to chord ratio for NACA airfoil' + ) + self.add_input('skin_thickness', val=0.002, units='m', desc='Skin panel thickness') + self.add_input( + 'twist_tail', + val=jnp.zeros(self.options['num_sections']), + units='deg', + desc='Twist distribution', + ) + + # Outputs + add_aviary_output( + self, Aircraft.VerticalTail.MASS, units='kg', desc='Total mass of the tail' + ) + + # File check + airfoil_file = self.options['airfoil_data_file'] + if airfoil_file is not None: + airfoil_file = get_path(airfoil_file) + airfoil_data = np.loadtxt(airfoil_file, skiprows=1) # Assume a header + + x_coords, y_coords = airfoil_data[:, 0], airfoil_data[:, 1] + + ( + self.camber, + self.camber_location, + self.max_thickness, + self.camber_line, + self.thickness, + ) = extract_airfoil_features(x_coords, y_coords) + + else: + NACA_digits = str(self.options['NACA_digits']) + # Parse the NACA airfoil type (4-digit) + self.camber = int(NACA_digits[0]) / 100.0 # Maximum camber + self.camber_location = int(NACA_digits[1]) / 10.0 # Location of max camber + self.max_thickness = int(NACA_digits[2:4]) / 100.0 # Max thickness + + def get_self_statics(self): + return ( + self.camber, + self.camber_location, + self.max_thickness, + self.camber_line, + self.thickness, + self.options['material'], + self.options['num_sections'], + self.options['NACA_digits'], + ) + + def compute_primal( + self, + aircraft__vertical_tail__span, + aircraft__vertical_tail__root_chord, + tip_chord_tail, + thickness_ratio, + skin_thickness, + twist_tail, + ): + material = self.options['material'] + density = materials.get_val(material, 'kg/m**3') + airfoil_file = self.options['airfoil_data_file'] + num_sections = self.options['num_sections'] + camber = self.camber + camber_location = self.camber_location + max_thickness = self.max_thickness + camber_line = self.camber_line + thickness = self.thickness + + # This is just so that the differentiation and unittest do not break. + aircraft__vertical_tail__mass = 0.0 * thickness_ratio + + # TODO unused?? + span_locations = jnp.linspace(0, aircraft__vertical_tail__span, num_sections) + + # Get x_points and dx for later + x_points, dx = precompute_airfoil_geometry(num_sections) + + # Thickness distribution + thickness_dist = airfoil_thickness(x_points, max_thickness) + + if airfoil_file is not None: + aircraft__vertical_tail__mass, _ = quadgk( + density * 2 * thickness * jnp.sqrt(1 + jnp.gradient(camber_line) ** 2), + [0, 1], + epsabs=1e-9, + epsrel=1e-9, + ) + else: + total_mass_first_part, _ = quadgk( + lambda x: density + * 2 + * jnp.atleast_1d(airfoil_thickness(x, max_thickness)) + * jnp.sqrt( + 1 + ((camber / camber_location**2) * (2 * camber_location - 2 * x)) ** 2 + ), + [0, camber_location], + epsabs=1e-9, + epsrel=1e-9, + ) + total_mass_second_part, _ = quadgk( + lambda x: density + * 2 + * jnp.atleast_1d(airfoil_thickness(x, max_thickness)) + * jnp.sqrt( + 1 + (camber / (1 - camber_location) ** 2 * (2 * camber_location - 2 * x)) ** 2 + ), + [camber_location, 1], + epsabs=1e-9, + epsrel=1e-9, + ) + + aircraft__vertical_tail__mass = total_mass_first_part + total_mass_second_part + + return aircraft__vertical_tail__mass diff --git a/aviary/subsystems/mass/simple_mass/test/__init__.py b/aviary/subsystems/mass/simple_mass/test/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/aviary/subsystems/mass/simple_mass/test/test_fuselage.py b/aviary/subsystems/mass/simple_mass/test/test_fuselage.py new file mode 100644 index 000000000..c27d4dacb --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/test/test_fuselage.py @@ -0,0 +1,49 @@ +import unittest + +import openmdao.api as om +from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal + +from aviary.subsystems.mass.simple_mass.fuselage import FuselageMass +from aviary.variable_info.variables import Aircraft + + +class FuselageMassTestCase(unittest.TestCase): + """Fuselage mass test case.""" + + def setUp(self): + self.prob = om.Problem() + self.prob.model.add_subsystem( + 'fuselage', + FuselageMass(), + promotes_inputs=['*'], + promotes_outputs=['*'], + ) + + self.prob.model.set_input_defaults(Aircraft.Fuselage.LENGTH, val=2.0, units='m') + + self.prob.model.set_input_defaults('base_diameter', val=0.5, units='m') + + self.prob.model.set_input_defaults('tip_diameter', val=0.3) + + self.prob.model.set_input_defaults('curvature', val=0.0, units='m') + + self.prob.model.set_input_defaults('y_offset', val=0.0, units='m') + + self.prob.model.set_input_defaults('z_offset', val=0.0, units='m') + + self.prob.setup(check=False, force_alloc_complex=True) + + def test_case(self): + self.prob.run_model() + + tol = 1e-3 + + assert_near_equal(self.prob[Aircraft.Fuselage.MASS], 373.849, tol) + + partial_data = self.prob.check_partials(out_stream=None, method='cs') + + assert_check_partials(partial_data, atol=1e-12, rtol=1e-12) + + +if __name__ == '__main__': + unittest.main() diff --git a/aviary/subsystems/mass/simple_mass/test/test_mass_summation.py b/aviary/subsystems/mass/simple_mass/test/test_mass_summation.py new file mode 100644 index 000000000..b5245e0cb --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/test/test_mass_summation.py @@ -0,0 +1,71 @@ +import unittest + +import openmdao.api as om +from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal + +from aviary.subsystems.mass.simple_mass.mass_summation import SimpleMassSummation, StructureMass +from aviary.variable_info.variables import Aircraft + + +class MassSummationTest(unittest.TestCase): + """Total mass summation test case.""" + + def test_case(self): + self.prob = om.Problem() + + self.prob.model.add_subsystem('tot', SimpleMassSummation(), promotes=['*']) + + self.prob.setup() + + self.prob.set_val(Aircraft.Fuselage.MASS, val=100.0) + self.prob.set_val(Aircraft.Wing.MASS, val=4.2) + self.prob.set_val(Aircraft.HorizontalTail.MASS, val=4.25) + self.prob.set_val(Aircraft.VerticalTail.MASS, val=4.5) + + self.prob.run_model() + + # om.n2(self.prob) + + tol = 1e-10 + + assert_near_equal(self.prob[Aircraft.Design.STRUCTURE_MASS], 112.95, tol) + + partial_data = self.prob.check_partials(out_stream=None, method='cs') + + assert_check_partials(partial_data) + + +class StructureMassTest(unittest.TestCase): + """Total structure summation mass test case.""" + + def setUp(self): + self.prob = om.Problem() + + self.prob.model.add_subsystem( + 'tot', + StructureMass(), + promotes_inputs=['*'], + promotes_outputs=['*'], + ) + + self.prob.setup(check=False, force_alloc_complex=True) + + self.prob.set_val(Aircraft.Fuselage.MASS, val=100.0) + self.prob.set_val(Aircraft.Wing.MASS, val=4.2) + self.prob.set_val(Aircraft.HorizontalTail.MASS, val=4.25) + self.prob.set_val(Aircraft.VerticalTail.MASS, val=4.5) + + def test_case(self): + self.prob.run_model() + + tol = 1e-10 + + assert_near_equal(self.prob[Aircraft.Design.STRUCTURE_MASS], 112.95, tol) + + partial_data = self.prob.check_partials(out_stream=None, method='cs') + + assert_check_partials(partial_data) + + +if __name__ == '__main__': + unittest.main() diff --git a/aviary/subsystems/mass/simple_mass/test/test_tail.py b/aviary/subsystems/mass/simple_mass/test/test_tail.py new file mode 100644 index 000000000..66272de98 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/test/test_tail.py @@ -0,0 +1,72 @@ +import unittest + +import numpy as np +import openmdao.api as om +from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal + +from aviary.subsystems.mass.simple_mass.tail import HorizontalTailMass, VerticalTailMass +from aviary.variable_info.variables import Aircraft + + +class TailMassTestCase(unittest.TestCase): + """Tail mass test case.""" + + def test_horizontal_tail(self): + prob = om.Problem() + prob.model.add_subsystem( + 'Tail', + HorizontalTailMass(), + promotes_inputs=['*'], + promotes_outputs=['*'], + ) + + prob.model.set_input_defaults(Aircraft.HorizontalTail.SPAN, val=1, units='m') + prob.model.set_input_defaults(Aircraft.HorizontalTail.ROOT_CHORD, val=1, units='m') + prob.model.set_input_defaults('tip_chord_tail', val=0.5, units='m') + prob.model.set_input_defaults('thickness_ratio', val=0.12) + prob.model.set_input_defaults('skin_thickness', val=0.002, units='m') + prob.model.set_input_defaults('twist_tail', val=np.zeros(10), units='deg') + + prob.setup(check=False, force_alloc_complex=True) + + prob.run_model() + + tol = 1e-4 + + assert_near_equal(prob[Aircraft.HorizontalTail.MASS], 10.6966719, tol) + + partial_data = prob.check_partials(out_stream=None, method='cs') + + assert_check_partials(partial_data, atol=1e-15, rtol=1e-15) + + def test_vertical_tail(self): + prob = om.Problem() + prob.model.add_subsystem( + 'Tail', + VerticalTailMass(), + promotes_inputs=['*'], + promotes_outputs=['*'], + ) + + prob.model.set_input_defaults(Aircraft.VerticalTail.SPAN, val=1, units='m') + prob.model.set_input_defaults(Aircraft.VerticalTail.ROOT_CHORD, val=1, units='m') + prob.model.set_input_defaults('tip_chord_tail', val=0.5, units='m') + prob.model.set_input_defaults('thickness_ratio', val=0.12) + prob.model.set_input_defaults('skin_thickness', val=0.002, units='m') + prob.model.set_input_defaults('twist_tail', val=np.zeros(10), units='deg') + + prob.setup(check=False, force_alloc_complex=True) + + prob.run_model() + + tol = 1e-4 + + assert_near_equal(prob[Aircraft.VerticalTail.MASS], 10.6966719, tol) + + partial_data = prob.check_partials(out_stream=None, method='cs') + + assert_check_partials(partial_data, atol=1e-15, rtol=1e-15) + + +if __name__ == '__main__': + unittest.main() diff --git a/aviary/subsystems/mass/simple_mass/test/test_wing.py b/aviary/subsystems/mass/simple_mass/test/test_wing.py new file mode 100644 index 000000000..45036ad83 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/test/test_wing.py @@ -0,0 +1,57 @@ +import unittest + +import jax.numpy as jnp +import numpy as np +import openmdao.api as om +from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal + +from aviary.subsystems.mass.simple_mass.wing import WingMass +from aviary.variable_info.variables import Aircraft + + +# @av.skipIfMissingDependencies(WingMassAndCOG) +class WingMassTestCase(unittest.TestCase): + """Wing mass test case.""" + + def setUp(self): + self.prob = om.Problem() + self.prob.model.add_subsystem( + 'wing_mass', + WingMass(), + promotes_inputs=['*'], + promotes_outputs=['*'], + ) + + self.prob.model.set_input_defaults(Aircraft.Wing.SPAN, val=1, units='m') + + self.prob.model.set_input_defaults(Aircraft.Wing.ROOT_CHORD, val=1, units='m') + + self.prob.model.set_input_defaults('tip_chord', val=0.5, units='m') + + self.prob.model.set_input_defaults('twist', val=jnp.zeros(10), units='deg') + + n_points = 10 # = num_sections + x = jnp.linspace(0, 1, n_points) + max_thickness_chord_ratio = 0.12 + thickness_dist = ( + 5 + * max_thickness_chord_ratio + * (0.2969 * np.sqrt(x) - 0.1260 * x - 0.3516 * x**2 + 0.2843 * x**3 - 0.1015 * x**4) + ) + + self.prob.model.set_input_defaults('thickness_dist', val=thickness_dist, units='m') + + self.prob.setup(check=False, force_alloc_complex=True) + + def test_case(self): + self.prob.run_model() + + tol = 1e-9 + assert_near_equal(self.prob[Aircraft.Wing.MASS], 10.6966719, tol) + + partial_data = self.prob.check_partials(out_stream=None, method='cs') + assert_check_partials(partial_data) + + +if __name__ == '__main__': + unittest.main() diff --git a/aviary/subsystems/mass/simple_mass/utils.py b/aviary/subsystems/mass/simple_mass/utils.py new file mode 100644 index 000000000..85ad74524 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/utils.py @@ -0,0 +1,58 @@ +import jax.numpy as jnp +import openmdao.jax as omj +from scipy.interpolate import CubicSpline + + +def precompute_airfoil_geometry(num_sections): + x_points = jnp.linspace(0, 1, num_sections) + dx = 1 / (num_sections - 1) + return x_points, dx + + +def airfoil_thickness(x, max_thickness): + return ( + 5 + * max_thickness + * (0.2969 * jnp.sqrt(x) - 0.1260 * x - 0.3516 * x**2 + 0.2843 * x**3 - 0.1015 * x**4) + ) + + +def airfoil_camber_line(x, camber, camber_location): + camber_location = omj.ks_max(camber_location, 1e-9) # Divide by zero check + return jnp.where( + x < camber_location, + (camber / camber_location**2) * (2 * camber_location * x - x**2), + (camber / (1 - camber_location) ** 2) + * ((1 - 2 * camber_location) + 2 * camber_location * x - x**2), + ) + + +def extract_airfoil_features(x_coords, y_coords): + """ + Extract camber, camber location, and max thickness from the given airfoil data. + This method assumes x_coords are normalized (ranging from 0 to 1). + """ + # Approximate the camber line and max thickness from the data + # Assume the camber line is the line of symmetry between the upper and lower surfaces + upper_surface = y_coords[: int(len(x_coords) // 2)] + lower_surface = y_coords[int(len(x_coords) // 2) :] + x_upper = x_coords[: int(len(x_coords) // 2)] + x_lower = x_coords[int(len(x_coords) // 2) :] + + upper_spline = CubicSpline(x_upper, upper_surface, bc_type='natural') + lower_spline = CubicSpline(x_lower, lower_surface, bc_type='natural') + + camber_line = (upper_spline(x_coords) + lower_spline(x_coords)) / 2.0 + + thickness = upper_spline(x_coords) - lower_spline(x_coords) + + max_thickness_index = omj.ks_max(thickness) + max_thickness_value = thickness[max_thickness_index] + + camber_slope = jnp.gradient(camber_line, x_coords) + camber_location_index = omj.ks_max(omj.smooth_abs(camber_slope)) + camber_location = x_coords[camber_location_index] + + camber = camber_line[camber_location_index] + + return camber, camber_location, max_thickness_value, camber_line, thickness diff --git a/aviary/subsystems/mass/simple_mass/wing.py b/aviary/subsystems/mass/simple_mass/wing.py new file mode 100644 index 000000000..ee44774e6 --- /dev/null +++ b/aviary/subsystems/mass/simple_mass/wing.py @@ -0,0 +1,168 @@ +from pathlib import Path + +import jax.numpy as jnp +import numpy as np +import openmdao.api as om + +from aviary.subsystems.mass.simple_mass.materials_database import materials +from aviary.subsystems.mass.simple_mass.utils import ( + airfoil_camber_line, + airfoil_thickness, + extract_airfoil_features, + precompute_airfoil_geometry, +) +from aviary.utils.functions import add_aviary_input, add_aviary_output, get_path +from aviary.utils.named_values import get_keys +from aviary.variable_info.variables import Aircraft + +try: + from quadax import quadgk +except ImportError: + raise ImportError( + "quadax package not found. You can install it by running 'pip install quadax'." + ) + + +class WingMass(om.JaxExplicitComponent): + def initialize(self): + self.options.declare('num_sections', types=int, default=10) + + self.options.declare( + 'NACA_digits', + default='2412', + types=(str, int), + desc='4 digit code for NACA airfoil of tail', + ) + + self.options.declare( + 'airfoil_data_file', + default=None, + types=(str, Path), + desc='File path for airfoil coordinates (overrides NACA digits)', + ) + + self.options.declare('num_sections', default=10, desc='Number of sections for enumeration') + + self.options.declare('material', default='Balsa', values=list(get_keys(materials))) + + self.camber = 0 + self.camber_location = 0 + self.max_thickness = 0 + self.camber_line = 0 + self.thickness = 0 + + def setup(self): + # Inputs + add_aviary_input(self, Aircraft.Wing.SPAN, units='m') # Full wingspan (adjustable) + + add_aviary_input(self, Aircraft.Wing.ROOT_CHORD, units='m') # Root chord length + + self.add_input('tip_chord', val=1.0, units='m') # Tip chord length -- no aviary input + + self.add_input( + 'twist', val=jnp.zeros(self.options['num_sections']), units='deg' + ) # Twist angles -- no aviary input + + self.add_input( + 'thickness_dist', + val=jnp.ones(self.options['num_sections']) * 0.1, + shape=(self.options['num_sections'],), + units='m', + ) # Thickness distribution of the wing (height) -- no aviary input + + # Outputs + add_aviary_output(self, Aircraft.Wing.MASS, units='kg') + + airfoil_data_file = self.options['airfoil_data_file'] + if airfoil_data_file is not None: + airfoil_data_file = get_path(airfoil_data_file) + airfoil_data = np.loadtxt(airfoil_data_file) + + x_coords, y_coords = airfoil_data[:, 0], airfoil_data[:, 1] + + ( + self.camber, + self.camber_location, + self.max_thickness, + self.thickness, + self.camber_line, + ) = extract_airfoil_features(x_coords, y_coords) + + else: + NACA_digits = str(self.options['NACA_digits']) + # Parse the NACA airfoil type (4-digit) + self.camber = int(NACA_digits[0]) / 100.0 # Maximum camber + self.camber_location = int(NACA_digits[1]) / 10.0 # Location of max camber + self.max_thickness = int(NACA_digits[2:4]) / 100.0 # Max thickness + + def get_self_statics(self): + return ( + self.camber, + self.camber_location, + self.max_thickness, + self.camber_line, + self.thickness, + self.options['material'], + self.options['num_sections'], + self.options['NACA_digits'], + ) + + def compute_primal( + self, aircraft__wing__span, aircraft__wing__root_chord, tip_chord, twist, thickness_dist + ): + material = self.options['material'] + density = materials.get_val(material, 'kg/m**3') + airfoil_data_file = self.options['airfoil_data_file'] + num_sections = self.options['num_sections'] + camber = self.camber + camber_location = self.camber_location + max_thickness = self.max_thickness + camber_line = self.camber_line + thickness = self.thickness + + # Get material density + density = materials.get_val(material, 'kg/m**3') + + # Wing spanwise distribution + jnp.linspace(0, aircraft__wing__span, num_sections) + + n_points = num_sections + jnp.linspace(0, 1, n_points) + 1 / (n_points - 1) + + if airfoil_data_file is not None: + aircraft__wing__mass, _ = quadgk( + density * 2 * thickness_dist * jnp.sqrt(1 + jnp.gradient(camber_line) ** 2), + [0, 1], + epsabs=1e-9, + epsrel=1e-9, + ) + else: + aircraft__wing__mass_first_part, _ = quadgk( + lambda x: density + * 2 + * jnp.atleast_1d(airfoil_thickness(x, max_thickness)) + * jnp.sqrt( + 1 + ((camber / camber_location**2) * (2 * camber_location - 2 * x)) ** 2 + ), + [0, camber_location], + epsabs=1e-9, + epsrel=1e-9, + ) + aircraft__wing__mass_second_part, _ = quadgk( + lambda x: density + * 2 + * jnp.atleast_1d(airfoil_thickness(x, max_thickness)) + * jnp.sqrt( + 1 + (camber / (1 - camber_location) ** 2 * (2 * camber_location - 2 * x)) ** 2 + ), + [camber_location, 1], + epsabs=1e-9, + epsrel=1e-9, + ) + + aircraft__wing__mass = ( + aircraft__wing__mass_first_part + aircraft__wing__mass_second_part + ) + + return aircraft__wing__mass diff --git a/pyproject.toml b/pyproject.toml index 34b98752f..a51b8abdd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,6 +12,7 @@ dependencies = [ "dymos>=1.14.0", "hvplot", "importlib_resources", + "jax", "matplotlib", "numpy<2", "openmdao>=3.37.0", @@ -19,6 +20,7 @@ dependencies = [ "panel>=1.0.0", "parameterized", "simupy", + "quadax" ] [project.optional-dependencies]