Simple Mass Component

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,

aviary/subsystems/mass/simple_mass/__init__.py

@@ -0,0 +1 @@
+

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

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()

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=['*'])

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

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')