Allow Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL to be in [1, 3]

aviary/subsystems/mass/flops_based/test/test_wing_detailed.py

@@ -6,6 +6,7 @@
 from openmdao.utils.testing_utils import use_tempdirs
 from parameterized import parameterized
 
+
 from aviary.subsystems.mass.flops_based.wing_detailed import (
     BWBDetailedWingBendingFact,
     DetailedWingBendingFact,
@@ -21,6 +22,7 @@
     print_case,
 )
 from aviary.variable_info.functions import setup_model_options
+from aviary.variable_info.options import get_option_defaults
 from aviary.variable_info.variables import Aircraft, Mission, Settings
 
 omit_cases = ['LargeSingleAisle2FLOPS']
@@ -34,7 +36,7 @@ class DetailedWingBendingTest(unittest.TestCase):
     def setUp(self):
         self.prob = om.Problem()
 
-    # Skip model that doesn't use detailed wing and BWB cases.
+    # Skip models that don't use detailed wing and skip BWB cases.
     @parameterized.expand(get_flops_case_names(omit=omit_cases), name_func=print_case)
     def test_case(self, case_name):
         prob = self.prob
@@ -420,6 +422,101 @@ def test_extreme_engine_loc(self):
         pod_inertia_factor = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
         assert_near_equal(pod_inertia_factor, 1.0, tolerance=1e-10)
 
+    def test_intensity_factor(self):
+        """
+        data taken from high_wing_single_aisle.csv
+        Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL = 3, 1, 1.5, 2.5
+        """
+        options = get_option_defaults()
+        options.set_val(Aircraft.Engine.NUM_ENGINES, val=[2], units='unitless')
+        options.set_val(Aircraft.Engine.NUM_WING_ENGINES, val=[2], units='unitless')
+        options.set_val(Aircraft.Propulsion.TOTAL_NUM_WING_ENGINES, val=2, units='unitless')
+        options.set_val(
+            Aircraft.Wing.INPUT_STATION_DISTRIBUTION, val=[0.0, 0.5, 1.0], units='unitless'
+        )
+        options.set_val(Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL, val=3, units='unitless')
+        options.set_val(Aircraft.Wing.NUM_INTEGRATION_STATIONS, val=46, units='unitless')
+
+        prob = self.prob
+        self.prob.model.add_subsystem(
+            'wing',
+            DetailedWingBendingFact(),
+            promotes_inputs=['*'],
+            promotes_outputs=['*'],
+        )
+
+        prob.model.set_input_defaults(
+            Aircraft.Wing.LOAD_PATH_SWEEP_DISTRIBUTION, val=[23, 22.0], units='deg'
+        )
+        prob.model.set_input_defaults(
+            Aircraft.Wing.THICKNESS_TO_CHORD_DISTRIBUTION, val=[0.1, 0.1, 0.1], units='unitless'
+        )
+        prob.model.set_input_defaults(
+            Aircraft.Wing.CHORD_PER_SEMISPAN_DISTRIBUTION, val=[0.13, 0.115, 0.06], units='unitless'
+        )
+        prob.model.set_input_defaults(Aircraft.Design.GROSS_MASS, val=150000.0, units='lbm')
+        prob.model.set_input_defaults(Aircraft.Wing.ASPECT_RATIO, val=9.40817168, units='unitless')
+        prob.model.set_input_defaults(
+            Aircraft.Wing.ASPECT_RATIO_REFERENCE, val=9.40817168, units='unitless'
+        )  # it was 0.01 in high_wing_single_aisle.csv
+        prob.model.set_input_defaults(
+            Aircraft.Wing.STRUT_BRACING_FACTOR, val=1.01, units='unitless'
+        )
+        prob.model.set_input_defaults(
+            Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR, val=0.01, units='unitless'
+        )
+        prob.model.set_input_defaults(Aircraft.Wing.THICKNESS_TO_CHORD, val=0.1, units='unitless')
+        prob.model.set_input_defaults(
+            Aircraft.Wing.THICKNESS_TO_CHORD_REFERENCE, val=0.01, units='unitless'
+        )
+
+        setup_model_options(prob, options)
+        prob.setup(check=False, force_alloc_complex=True)
+
+        prob.set_val(Aircraft.Engine.POD_MASS, np.array([5587.40886136]), units='lbm')
+        prob.set_val(Aircraft.Engine.WING_LOCATIONS, [0.22])
+
+        prob.run_model()
+
+        bending_mat_factor = prob.get_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR)
+        assert_near_equal(bending_mat_factor, 2.9367664396, tolerance=1e-9)
+
+        pod_inertia_factor = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
+        assert_near_equal(pod_inertia_factor, 0.9772998541, tolerance=1e-9)
+
+        options.set_val(Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL, val=1, units='unitless')
+        setup_model_options(prob, options)
+        prob.setup(check=False, force_alloc_complex=True)
+        prob.run_model()
+
+        bending_mat_factor = prob.get_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR)
+        assert_near_equal(bending_mat_factor, 1.51247369, tolerance=1e-9)
+
+        pod_inertia_factor = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
+        assert_near_equal(pod_inertia_factor, 1.0, tolerance=1e-9)
+
+        options.set_val(Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL, val=1.5, units='unitless')
+        setup_model_options(prob, options)
+        prob.setup(check=False, force_alloc_complex=True)
+        prob.run_model()
+
+        bending_mat_factor = prob.get_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR)
+        assert_near_equal(bending_mat_factor, 1.94203494, tolerance=1e-9)
+
+        pod_inertia_factor = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
+        assert_near_equal(pod_inertia_factor, 1.0, tolerance=1e-9)
+
+        options.set_val(Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL, val=2.5, units='unitless')
+        setup_model_options(prob, options)
+        prob.setup(check=False, force_alloc_complex=True)
+        prob.run_model()
+
+        bending_mat_factor = prob.get_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR)
+        assert_near_equal(bending_mat_factor, 2.61652282, tolerance=1e-9)
+
+        pod_inertia_factor = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
+        assert_near_equal(pod_inertia_factor, 1.0, tolerance=1e-9)
+
     def test_IO(self):
         assert_match_varnames(self.prob.model)
 
@@ -607,7 +704,8 @@ def test_case1(self):
         pod_inertia_expected = 1.0
         assert_near_equal(BENDING_MATERIAL_FACTOR, BENDING_MATERIAL_FACTOR_expected, tolerance=1e-9)
         assert_near_equal(prob.get_val('calculated_wing_area'), 5399.4057051, tolerance=1e-9)
-        # current BWB data set does not check the following
+        # For BWB with non wing engines, inertia_factor_prod = 1 always.
+        # Need a different dataset with wing engines to check this one
         assert_near_equal(pod_inertia, pod_inertia_expected, tolerance=1e-9)
 
     def test_case2(self):
@@ -693,12 +791,11 @@ def test_case2(self):
         pod_inertia_expected = 1.0
         assert_near_equal(BENDING_MATERIAL_FACTOR, BENDING_MATERIAL_FACTOR_expected, tolerance=1e-9)
         assert_near_equal(prob.get_val('calculated_wing_area'), 4151.88659141, tolerance=1e-9)
-        # current BWB data set does not check the following
         assert_near_equal(pod_inertia, pod_inertia_expected, tolerance=1e-9)
 
 
 if __name__ == '__main__':
     unittest.main()
     # test = DetailedWingBendingTest()
     # test.setUp()
-    # test.test_extreme_engine_loc()
+    # test.test_intensity_factor()

aviary/subsystems/mass/flops_based/wing_detailed.py

@@ -5,7 +5,45 @@
 
 from aviary.variable_info.enums import Verbosity
 from aviary.variable_info.functions import add_aviary_input, add_aviary_option, add_aviary_output
-from aviary.variable_info.variables import Aircraft, Mission, Settings
+from aviary.variable_info.variables import Aircraft, Settings
+
+
+def load_intensity_by_factor(load_dist_factor, intn_stations):
+    """
+    Calculate load intensity at wing integration stations for the given load_dist_factor.
+
+    Parameters:
+    load_dist_factor (float): 0 or 1 <= load_dist_factor <= 3.
+        1.0 : triangular distribution
+        2.0 : elliptical distribution (default)
+        3.0 : rectangular distribution
+        1.0-2.0 : blend of triangular and elliptical
+        2.0-3.0 : blend of elliptical and rectangular
+    intn_stations (array): integration stations
+
+    Note: In FLOPS, load_dist_factor can be 0 in which case load intensities are computed
+          based on input pressure distribution.
+    """
+    if load_dist_factor == 1:
+        load_intensity = 1.0 - intn_stations
+    elif load_dist_factor == 2:
+        load_intensity = np.sqrt(1.0 - intn_stations**2)
+    elif load_dist_factor == 3:
+        load_intensity = np.ones(len(intn_stations))
+    elif 1 < load_dist_factor < 2:
+        load_intensity = (load_dist_factor - 1.0) * np.sqrt(1.0 - intn_stations**2) + (
+            2.0 - load_dist_factor
+        ) * (1.0 - intn_stations)
+    elif 2 < load_dist_factor < 3:
+        load_intensity = (3.0 - load_dist_factor) * np.sqrt(1.0 - intn_stations**2) + (
+            load_dist_factor - 2.0
+        ) * np.ones(len(intn_stations))
+    else:
+        raise om.AnalysisError(
+            f'{load_dist_factor} is not a valid value for '
+            f'{Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL}, it must be between 1 and 3 (inclusive).'
+        )
+    return load_intensity
 
 
 class DetailedWingBendingFact(om.ExplicitComponent):
@@ -85,15 +123,6 @@ def compute(self, inputs, outputs):
         num_wing_engines = self.options[Aircraft.Engine.NUM_WING_ENGINES]
         num_engine_type = len(num_wing_engines)
 
-        # TODO: Support all options for this parameter.
-        # 0.0 : input distribution
-        # 1.0 : triangular distribution
-        # 2.0 : elliptical distribution (default)
-        # 3.0 : rectangular distribution
-        # 1.0-2.0 : blend of triangular and elliptical
-        # 2.0-3.0 : blend of elliptical and rectangular
-        load_distribution_factor = self.options[Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL]
-
         load_path_sweep = inputs[Aircraft.Wing.LOAD_PATH_SWEEP_DISTRIBUTION]
         thickness_to_chord = inputs[Aircraft.Wing.THICKNESS_TO_CHORD_DISTRIBUTION]
         chord = inputs[Aircraft.Wing.CHORD_PER_SEMISPAN_DISTRIBUTION]
@@ -138,18 +167,8 @@ def compute(self, inputs, outputs):
             (dy[1:] + 2.0 * integration_stations[1:-1]) * dy[1:] * sweep_int_stations[1:-1]
         )
 
-        # TODO: add all load_distribution_factor options
-        if load_distribution_factor == 1:
-            load_intensity = 1.0 - integration_stations
-        elif load_distribution_factor == 2:
-            load_intensity = np.sqrt(1.0 - integration_stations**2)
-        elif load_distribution_factor == 3:
-            load_intensity = np.ones(num_integration_stations + 1)
-        else:
-            raise om.AnalysisError(
-                f'{load_distribution_factor} is not a valid value for '
-                f'{Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL}, it must be "1", "2", or "3".'
-            )
+        load_distribution_factor = self.options[Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL]
+        load_intensity = load_intensity_by_factor(load_distribution_factor, integration_stations)
 
         chord_interp = InterpND(
             method='slinear', points=(inp_stations), x_interp=integration_stations
@@ -341,8 +360,6 @@ def setup_partials(self):
     def compute(self, inputs, outputs):
         verbosity = self.options[Settings.VERBOSITY]
         num_integration_stations = self.options[Aircraft.Wing.NUM_INTEGRATION_STATIONS]
-        num_wing_engines = self.options[Aircraft.Engine.NUM_WING_ENGINES]
-        num_engine_type = len(num_wing_engines)
         width = inputs[Aircraft.Fuselage.MAX_WIDTH][0]
         wingspan = inputs[Aircraft.Wing.SPAN][0]
         rate_span = (wingspan - width) / wingspan
@@ -370,15 +387,6 @@ def compute(self, inputs, outputs):
         # For BWB, always start from inp_stations_mod[1], not inp_stations_mod[0]
         inp_stations_mod = inp_stations_mod[1:]
 
-        # TODO: Support all options for this parameter.
-        # 0.0 : input distribution
-        # 1.0 : triangular distribution
-        # 2.0 : elliptical distribution (default)
-        # 3.0 : rectangular distribution
-        # 1.0-2.0 : blend of triangular and elliptical
-        # 2.0-3.0 : blend of elliptical and rectangular
-        load_distribution_factor = self.options[Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL]
-
         load_path_sweep = inputs['BWB_LOAD_PATH_SWEEP_DISTRIBUTION']
         load_path_sweep_mod = np.array(load_path_sweep[1:])
 
@@ -439,18 +447,8 @@ def compute(self, inputs, outputs):
             (dy[0:] + 2.0 * integration_stations[0:-1]) * dy[0:] * sweep_int_stations[0:-1]
         )
 
-        # TODO: add all load_distribution_factor options
-        if load_distribution_factor == 1:
-            load_intensity = 1.0 - integration_stations
-        elif load_distribution_factor == 2:
-            load_intensity = np.sqrt(1.0 - integration_stations**2)
-        elif load_distribution_factor == 3:
-            load_intensity = np.ones(num_integration_stations + 1)
-        else:
-            raise om.AnalysisError(
-                f'{load_distribution_factor} is not a valid value for '
-                f'{Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL}, it must be "1", "2", or "3".'
-            )
+        load_distribution_factor = self.options[Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL]
+        load_intensity = load_intensity_by_factor(load_distribution_factor, integration_stations)
 
         chord_interp = InterpND(
             method='slinear', points=(inp_stations_mod), x_interp=integration_stations
@@ -523,13 +521,14 @@ def compute(self, inputs, outputs):
 
         outputs[Aircraft.Wing.BENDING_MATERIAL_FACTOR] = bt
 
+        num_wing_engines = self.options[Aircraft.Engine.NUM_WING_ENGINES]
+        num_engine_type = len(num_wing_engines)
         engine_locations = inputs[Aircraft.Engine.WING_LOCATIONS]
         gross_mass = inputs[Aircraft.Design.GROSS_MASS]
         # NOTE pod mass assumed the same for wing/non-wing mounted engines, only using
         #      wing mounted pods here
         pod_mass = inputs[Aircraft.Engine.POD_MASS]
         if np.sum(num_wing_engines) > 0:
-            # TODO: the rest is not checked.
             inertia_factor = np.zeros(num_engine_type, dtype=chord.dtype)
             eel = np.zeros(len(dy) + 1, dtype=chord.dtype)
 

aviary/variable_info/variable_meta_data.py

@@ -5690,7 +5690,7 @@
         'FLOPS': 'WTIN.PDIST',  # ['&DEFINE.WTIN.PDIST', 'WDEF.PDIST'],
     },
     units='unitless',
-    desc='controls spatial distribution of integration stations for detailed wing',
+    desc='controls spatial distribution of integration stations for detailed wing, in [1, 3]',
     default_value=2.0,
     option=True,
 )