diff --git a/aviary/docs/examples/modified_aircraft.csv b/aviary/docs/examples/modified_aircraft.csv
index d354187a8..43ac7ba2f 100644
--- a/aviary/docs/examples/modified_aircraft.csv
+++ b/aviary/docs/examples/modified_aircraft.csv
@@ -114,7 +114,7 @@ aircraft:wing:aeroelastic_tailoring_factor,0.0,unitless
 aircraft:wing:airfoil_technology,1.92669766647637,unitless
 aircraft:wing:area,1370.0,ft**2
 aircraft:wing:aspect_ratio,11.02091,unitless
-aircraft:wing:bending_mass_scaler,1.0,unitless
+aircraft:wing:bending_material_mass_scaler,1.0,unitless
 aircraft:wing:chord_per_semispan,0.31,0.23,0.084,unitless
 aircraft:wing:composite_fraction,0.2,unitless
 aircraft:wing:control_surface_area,137,ft**2
diff --git a/aviary/interface/methods_for_level2.py b/aviary/interface/methods_for_level2.py
index dc66c9e13..6c7df771f 100644
--- a/aviary/interface/methods_for_level2.py
+++ b/aviary/interface/methods_for_level2.py
@@ -334,6 +334,7 @@ def load_inputs(
         self.mass_method = mass_method = aviary_inputs.get_val(Settings.MASS_METHOD)
 
         if mission_method is TWO_DEGREES_OF_FREEDOM or mass_method is GASP:
+            # TODO this should be a preprocessor step if it is required here
             aviary_inputs = update_GASP_options(aviary_inputs)
         initialization_guesses = initialization_guessing(
             aviary_inputs, initialization_guesses, engine_builders
@@ -2799,8 +2800,8 @@ def _add_guesses(self, phase_name, phase, guesses, setvalprob, parent_prefix):
                 # are not consistent
                 if initial_bounds[1] != duration_bounds[1]:
                     raise ValueError(
-                        f"Initial and duration bounds for {
-                            phase_name} are not consistent."
+                        f"Initial and duration bounds for {phase_name} "
+                        "are not consistent."
                     )
                 guesses["time"] = (
                     [np.mean(initial_bounds[0]), np.mean(duration_bounds[0])],
@@ -2867,8 +2868,8 @@ def _add_guesses(self, phase_name, phase, guesses, setvalprob, parent_prefix):
                 else:
                     # raise error if the guess key is not recognized
                     raise ValueError(
-                        f"Initial guess key {guess_key} in {
-                            phase_name} is not recognized."
+                        f"Initial guess key {guess_key} in {phase_name} is not "
+                        "recognized."
                     )
 
         if self.mission_method is SOLVED_2DOF:
diff --git a/aviary/models/N3CC/N3CC_data.py b/aviary/models/N3CC/N3CC_data.py
index ae935688d..8c043da8d 100644
--- a/aviary/models/N3CC/N3CC_data.py
+++ b/aviary/models/N3CC/N3CC_data.py
@@ -247,7 +247,7 @@
 inputs.set_val(Aircraft.Wing.AREA, 1220.0, 'ft**2')
 inputs.set_val(Aircraft.Wing.ASPECT_RATIO, 11.5587605382765)
 inputs.set_val(Aircraft.Wing.ASPECT_RATIO_REF, 11.5587605382765)
-inputs.set_val(Aircraft.Wing.BENDING_MASS_SCALER, 1.0)
+inputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER, 1.0)
 inputs.set_val(
     Aircraft.Wing.CHORD_PER_SEMISPAN_DIST,
     np.array(
@@ -432,8 +432,8 @@
 outputs.set_val(Aircraft.VerticalTail.FINENESS, 0.1000)
 outputs.set_val(Aircraft.VerticalTail.MASS, 1175.0, 'lbm')
 
-outputs.set_val(Aircraft.Wing.BENDING_FACTOR, 11.9602)
-outputs.set_val(Aircraft.Wing.BENDING_MASS, 5410.5, 'lbm')
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR, 11.9602)
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS, 5410.5, 'lbm')
 outputs.set_val(Aircraft.Wing.CHARACTERISTIC_LENGTH, 10.27, 'ft')
 outputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA, 0.333 * 1220, 'ft**2')
 outputs.set_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR, 0.960516)
diff --git a/aviary/models/N3CC/N3CC_generic_low_speed_polars_FLOPSinp.csv b/aviary/models/N3CC/N3CC_generic_low_speed_polars_FLOPSinp.csv
index a222359a8..17d699c42 100644
--- a/aviary/models/N3CC/N3CC_generic_low_speed_polars_FLOPSinp.csv
+++ b/aviary/models/N3CC/N3CC_generic_low_speed_polars_FLOPSinp.csv
@@ -116,7 +116,7 @@ aircraft:wing:airfoil_technology,1.6,unitless
 aircraft:wing:area,1220,1,0,0,0,0,ft**2
 aircraft:wing:aspect_ratio,11.5587605382765,1,0,0,0,0,unitless
 aircraft:wing:aspect_ratio_reference,11.5587605382765,unitless
-aircraft:wing:bending_mass_scaler,1,unitless
+aircraft:wing:BENDING_MATERIAL_MASS_SCALER,1,unitless
 aircraft:wing:bwb_aft_body_mass_scaler,1,unitless
 aircraft:wing:chord_per_semispan,0.273522534166506,0.204274849507037,0.0888152947868224,0.0725353313595661,unitless
 aircraft:wing:composite_fraction,0.33333,unitless
diff --git a/aviary/models/large_single_aisle_1/large_single_aisle_1_FLOPS_data.py b/aviary/models/large_single_aisle_1/large_single_aisle_1_FLOPS_data.py
index 9c74cfbd5..5ec622a3f 100644
--- a/aviary/models/large_single_aisle_1/large_single_aisle_1_FLOPS_data.py
+++ b/aviary/models/large_single_aisle_1/large_single_aisle_1_FLOPS_data.py
@@ -194,7 +194,7 @@
 inputs.set_val(Aircraft.Wing.AIRFOIL_TECHNOLOGY, 1.92669766647637)
 inputs.set_val(Aircraft.Wing.AREA, 1370.0, 'ft**2')
 inputs.set_val(Aircraft.Wing.ASPECT_RATIO, 11.22091)
-inputs.set_val(Aircraft.Wing.BENDING_MASS_SCALER, 1.0)
+inputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER, 1.0)
 inputs.set_val(Aircraft.Wing.CHORD_PER_SEMISPAN_DIST, np.array([0.31, 0.23, 0.084]))
 inputs.set_val(Aircraft.Wing.COMPOSITE_FRACTION, 0.2)
 inputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA, 137, 'ft**2')
@@ -369,8 +369,8 @@
 outputs.set_val(Aircraft.VerticalTail.FINENESS, 0.1195)
 outputs.set_val(Aircraft.VerticalTail.MASS, 1221.8, 'lbm')
 
-outputs.set_val(Aircraft.Wing.BENDING_FACTOR, 11.5918)
-outputs.set_val(Aircraft.Wing.BENDING_MASS, 8184.8, 'lbm')
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR, 11.5918)
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS, 8184.8, 'lbm')
 outputs.set_val(Aircraft.Wing.CHARACTERISTIC_LENGTH, 10.49, 'ft')
 # Not in FLOPS output; calculated from inputs.
 outputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA, 137, 'ft**2')
diff --git a/aviary/models/large_single_aisle_2/large_single_aisle_2_FLOPS_data.py b/aviary/models/large_single_aisle_2/large_single_aisle_2_FLOPS_data.py
index 3787f9ddc..889dbc855 100644
--- a/aviary/models/large_single_aisle_2/large_single_aisle_2_FLOPS_data.py
+++ b/aviary/models/large_single_aisle_2/large_single_aisle_2_FLOPS_data.py
@@ -217,7 +217,7 @@
 inputs.set_val(Aircraft.Wing.AIRFOIL_TECHNOLOGY, 1.87)
 inputs.set_val(Aircraft.Wing.AREA, 1341.0, 'ft**2')
 inputs.set_val(Aircraft.Wing.ASPECT_RATIO, 9.45)
-inputs.set_val(Aircraft.Wing.BENDING_MASS_SCALER, 1.0)
+inputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER, 1.0)
 inputs.set_val(Aircraft.Wing.COMPOSITE_FRACTION, 0.0)
 inputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA_RATIO, 0.333)
 inputs.set_val(Aircraft.Wing.GLOVE_AND_BAT, 0.0, 'ft**2')
@@ -397,8 +397,8 @@
 outputs.set_val(Aircraft.VerticalTail.FINENESS, 0.1375)
 outputs.set_val(Aircraft.VerticalTail.MASS, 1035.6, 'lbm')
 
-outputs.set_val(Aircraft.Wing.BENDING_FACTOR, 8.8294)
-outputs.set_val(Aircraft.Wing.BENDING_MASS, 6016.9, 'lbm')
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR, 8.8294)
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS, 6016.9, 'lbm')
 outputs.set_val(Aircraft.Wing.CHARACTERISTIC_LENGTH, 11.91, 'ft')
 outputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA, 0.333 * 1341.0, 'ft**2')
 outputs.set_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR, 0.940000)
diff --git a/aviary/models/large_single_aisle_2/large_single_aisle_2_altwt_FLOPS_data.py b/aviary/models/large_single_aisle_2/large_single_aisle_2_altwt_FLOPS_data.py
index 454b605f4..322332a86 100644
--- a/aviary/models/large_single_aisle_2/large_single_aisle_2_altwt_FLOPS_data.py
+++ b/aviary/models/large_single_aisle_2/large_single_aisle_2_altwt_FLOPS_data.py
@@ -204,7 +204,7 @@
 inputs.set_val(Aircraft.Wing.AIRFOIL_TECHNOLOGY, 1.87)
 inputs.set_val(Aircraft.Wing.AREA, 1341.0, 'ft**2')
 inputs.set_val(Aircraft.Wing.ASPECT_RATIO, 9.45)
-inputs.set_val(Aircraft.Wing.BENDING_MASS_SCALER, 1.0)
+inputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER, 1.0)
 inputs.set_val(Aircraft.Wing.COMPOSITE_FRACTION, 0.0)
 inputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA_RATIO, 0.333)
 inputs.set_val(Aircraft.Wing.GLOVE_AND_BAT, 0.0, 'ft**2')
@@ -375,8 +375,8 @@
 outputs.set_val(Aircraft.VerticalTail.FINENESS, 0.1375)
 outputs.set_val(Aircraft.VerticalTail.MASS, 1707., 'lbm')
 
-outputs.set_val(Aircraft.Wing.BENDING_FACTOR, 8.8294)
-outputs.set_val(Aircraft.Wing.BENDING_MASS, 6016.9, 'lbm')
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR, 8.8294)
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS, 6016.9, 'lbm')
 outputs.set_val(Aircraft.Wing.CHARACTERISTIC_LENGTH, 11.91, 'ft')
 outputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA, 0.333 * 1341.0, 'ft**2')
 outputs.set_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR, 0.940000)
diff --git a/aviary/models/large_single_aisle_2/large_single_aisle_2_detailwing_FLOPS_data.py b/aviary/models/large_single_aisle_2/large_single_aisle_2_detailwing_FLOPS_data.py
index 63721e4a4..7711d695f 100644
--- a/aviary/models/large_single_aisle_2/large_single_aisle_2_detailwing_FLOPS_data.py
+++ b/aviary/models/large_single_aisle_2/large_single_aisle_2_detailwing_FLOPS_data.py
@@ -204,7 +204,7 @@
 inputs.set_val(Aircraft.Wing.AIRFOIL_TECHNOLOGY, 1.87)
 inputs.set_val(Aircraft.Wing.AREA, 1341.0, 'ft**2')
 inputs.set_val(Aircraft.Wing.ASPECT_RATIO, 9.42519)
-inputs.set_val(Aircraft.Wing.BENDING_MASS_SCALER, 1.0)
+inputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER, 1.0)
 inputs.set_val(
     Aircraft.Wing.CHORD_PER_SEMISPAN_DIST, np.array([0.4441, 0.2313, 0.0729]))
 inputs.set_val(Aircraft.Wing.COMPOSITE_FRACTION, 0.0)
@@ -394,8 +394,8 @@
 outputs.set_val(Aircraft.VerticalTail.FINENESS, 0.1375)
 outputs.set_val(Aircraft.VerticalTail.MASS, 1035.6, 'lbm')
 
-outputs.set_val(Aircraft.Wing.BENDING_FACTOR, 9.0236)
-outputs.set_val(Aircraft.Wing.BENDING_MASS, 6276.3, 'lbm')
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR, 9.0236)
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS, 6276.3, 'lbm')
 outputs.set_val(Aircraft.Wing.CHARACTERISTIC_LENGTH, 11.91, 'ft')
 outputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA, 0.333 * 1341.0, 'ft**2')
 outputs.set_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR, 0.959104)
diff --git a/aviary/models/multi_engine_single_aisle/multi_engine_single_aisle_data.py b/aviary/models/multi_engine_single_aisle/multi_engine_single_aisle_data.py
index 0ff5b633c..f5773cf40 100644
--- a/aviary/models/multi_engine_single_aisle/multi_engine_single_aisle_data.py
+++ b/aviary/models/multi_engine_single_aisle/multi_engine_single_aisle_data.py
@@ -232,7 +232,7 @@
 inputs.set_val(Aircraft.Wing.AIRFOIL_TECHNOLOGY, 1.87)
 inputs.set_val(Aircraft.Wing.AREA, 1341.0, 'ft**2')
 inputs.set_val(Aircraft.Wing.ASPECT_RATIO, 9.42519)
-inputs.set_val(Aircraft.Wing.BENDING_MASS_SCALER, 1.0)
+inputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER, 1.0)
 inputs.set_val(
     Aircraft.Wing.CHORD_PER_SEMISPAN_DIST, np.array([0.4441, 0.2313, 0.0729]))
 inputs.set_val(Aircraft.Wing.COMPOSITE_FRACTION, 0.0)
@@ -410,8 +410,8 @@
 outputs.set_val(Aircraft.VerticalTail.FINENESS, 0.1195)
 outputs.set_val(Aircraft.VerticalTail.MASS, 1221.8, 'lbm')
 
-outputs.set_val(Aircraft.Wing.BENDING_FACTOR, 11.5918)
-outputs.set_val(Aircraft.Wing.BENDING_MASS, 8184.8, 'lbm')
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR, 11.5918)
+outputs.set_val(Aircraft.Wing.BENDING_MATERIAL_MASS, 8184.8, 'lbm')
 outputs.set_val(Aircraft.Wing.CHARACTERISTIC_LENGTH, 10.49, 'ft')
 # Not in FLOPS output; calculated from inputs.
 outputs.set_val(Aircraft.Wing.CONTROL_SURFACE_AREA, 137, 'ft**2')
diff --git a/aviary/models/test_aircraft/aircraft_for_bench_FwFm.csv b/aviary/models/test_aircraft/aircraft_for_bench_FwFm.csv
index bd9fa3412..29b7a41ec 100644
--- a/aviary/models/test_aircraft/aircraft_for_bench_FwFm.csv
+++ b/aviary/models/test_aircraft/aircraft_for_bench_FwFm.csv
@@ -114,7 +114,7 @@ aircraft:wing:aeroelastic_tailoring_factor,0.0,unitless
 aircraft:wing:airfoil_technology,1.92669766647637,unitless
 aircraft:wing:area,1370.0,ft**2
 aircraft:wing:aspect_ratio,11.22091,unitless
-aircraft:wing:bending_mass_scaler,1.0,unitless
+aircraft:wing:BENDING_MATERIAL_MASS_SCALER,1.0,unitless
 aircraft:wing:chord_per_semispan,0.31,0.23,0.084,unitless
 aircraft:wing:composite_fraction,0.2,unitless
 aircraft:wing:control_surface_area,137,ft**2
diff --git a/aviary/models/test_aircraft/aircraft_for_bench_FwFm_with_electric.csv b/aviary/models/test_aircraft/aircraft_for_bench_FwFm_with_electric.csv
index 7dcd99288..4eee8d79b 100644
--- a/aviary/models/test_aircraft/aircraft_for_bench_FwFm_with_electric.csv
+++ b/aviary/models/test_aircraft/aircraft_for_bench_FwFm_with_electric.csv
@@ -114,7 +114,7 @@ aircraft:wing:aeroelastic_tailoring_factor,0.0,unitless
 aircraft:wing:airfoil_technology,1.92669766647637,unitless
 aircraft:wing:area,1370.0,ft**2
 aircraft:wing:aspect_ratio,11.22091,unitless
-aircraft:wing:bending_mass_scaler,1.0,unitless
+aircraft:wing:BENDING_MATERIAL_MASS_SCALER,1.0,unitless
 aircraft:wing:chord_per_semispan,0.31,0.23,0.084,unitless
 aircraft:wing:composite_fraction,0.2,unitless
 aircraft:wing:control_surface_area,137,ft**2
diff --git a/aviary/models/test_aircraft/aircraft_for_bench_FwGm.csv b/aviary/models/test_aircraft/aircraft_for_bench_FwGm.csv
index e209ed2ca..8d97f49aa 100644
--- a/aviary/models/test_aircraft/aircraft_for_bench_FwGm.csv
+++ b/aviary/models/test_aircraft/aircraft_for_bench_FwGm.csv
@@ -219,7 +219,7 @@ aircraft:vertical_tail:mass_scaler,1.0,unitless
 aircraft:vertical_tail:num_tails,1,unitless
 aircraft:wing:aeroelastic_tailoring_factor,0.0,unitless
 aircraft:wing:airfoil_technology,1.92669766647637,unitless
-aircraft:wing:bending_mass_scaler,1.0,unitless
+aircraft:wing:BENDING_MATERIAL_MASS_SCALER,1.0,unitless
 aircraft:wing:chord_per_semispan,0.31,0.23,0.084,unitless
 aircraft:wing:composite_fraction,0.2,unitless
 aircraft:wing:control_surface_area,137,ft**2
diff --git a/aviary/models/test_aircraft/aircraft_for_bench_GwFm.csv b/aviary/models/test_aircraft/aircraft_for_bench_GwFm.csv
index 092cfcf3e..4b39a7868 100644
--- a/aviary/models/test_aircraft/aircraft_for_bench_GwFm.csv
+++ b/aviary/models/test_aircraft/aircraft_for_bench_GwFm.csv
@@ -204,7 +204,7 @@ aircraft:wing:airfoil_technology,1.92669766647637,unitless
 aircraft:wing:area,1370.0,ft**2
 aircraft:wing:aspect_ratio,10.13,unitless
 aircraft:wing:average_chord,12.615,ft
-aircraft:wing:bending_mass_scaler,1.0,unitless
+aircraft:wing:BENDING_MATERIAL_MASS_SCALER,1.0,unitless
 aircraft:wing:center_distance,0.463,unitless
 aircraft:wing:choose_fold_location,True,unitless
 aircraft:wing:chord_per_semispan,0.31,0.23,0.084,unitless
diff --git a/aviary/models/test_aircraft/aircraft_for_bench_solved2dof.csv b/aviary/models/test_aircraft/aircraft_for_bench_solved2dof.csv
index b6848994f..3aebc639f 100644
--- a/aviary/models/test_aircraft/aircraft_for_bench_solved2dof.csv
+++ b/aviary/models/test_aircraft/aircraft_for_bench_solved2dof.csv
@@ -114,7 +114,7 @@ aircraft:wing:aeroelastic_tailoring_factor,0.0,unitless
 aircraft:wing:airfoil_technology,1.92669766647637,unitless
 aircraft:wing:area,1370.0,ft**2
 aircraft:wing:aspect_ratio,11.22091,unitless
-aircraft:wing:bending_mass_scaler,1.0,unitless
+aircraft:wing:BENDING_MATERIAL_MASS_SCALER,1.0,unitless
 aircraft:wing:chord_per_semispan,0.31,0.23,0.084,unitless
 aircraft:wing:composite_fraction,0.2,unitless
 aircraft:wing:control_surface_area,137,ft**2
diff --git a/aviary/subsystems/aerodynamics/flops_based/test/data/high_wing_single_aisle.csv b/aviary/subsystems/aerodynamics/flops_based/test/data/high_wing_single_aisle.csv
index ebde1eae7..d798dfe0b 100644
--- a/aviary/subsystems/aerodynamics/flops_based/test/data/high_wing_single_aisle.csv
+++ b/aviary/subsystems/aerodynamics/flops_based/test/data/high_wing_single_aisle.csv
@@ -98,7 +98,7 @@ aircraft:wing:aeroelastic_tailoring_factor,0.01,unitless   #check
 aircraft:wing:airfoil_technology,1.01,unitless   #check
 aircraft:wing:area,1480,ft**2
 aircraft:wing:aspect_ratio_reference,0.01,unitless   #check
-aircraft:wing:bending_mass_scaler,1.01,unitless   #check
+aircraft:wing:BENDING_MATERIAL_MASS_SCALER,1.01,unitless   #check
 aircraft:wing:chord_per_semispan,0.13,0.115,0.06,unitless
 aircraft:wing:composite_fraction,0.01,unitless   #check
 aircraft:wing:dihedral,-1.0,deg
diff --git a/aviary/subsystems/mass/flops_based/landing_gear.py b/aviary/subsystems/mass/flops_based/landing_gear.py
index db87a4aee..f053a2586 100644
--- a/aviary/subsystems/mass/flops_based/landing_gear.py
+++ b/aviary/subsystems/mass/flops_based/landing_gear.py
@@ -294,8 +294,12 @@ def setup(self):
         add_aviary_input(self, Aircraft.Fuselage.MAX_WIDTH, val=0.0)
         add_aviary_input(self, Aircraft.Nacelle.AVG_DIAMETER,
                          val=np.zeros(num_engine_type))
-        add_aviary_input(self, Aircraft.Engine.WING_LOCATIONS,
-                         val=np.zeros((num_engine_type, int(num_wing_engines[0]/2))))
+        if any(num_wing_engines) > 0:
+            add_aviary_input(self, Aircraft.Engine.WING_LOCATIONS, val=np.zeros(
+                (num_engine_type, int(num_wing_engines[0] / 2))))
+        else:
+            add_aviary_input(self, Aircraft.Engine.WING_LOCATIONS,
+                             val=[[0.0]])
         add_aviary_input(self, Aircraft.Wing.DIHEDRAL, val=0.0)
         add_aviary_input(self, Aircraft.Wing.SPAN, val=0.0)
 
@@ -309,7 +313,6 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         # TODO temp using first engine, heterogeneous engines not supported
         num_eng = options.get_val(Aircraft.Engine.NUM_ENGINES)[0]
         num_wing_eng = options.get_val(Aircraft.Engine.NUM_WING_ENGINES)[0]
-
         y_eng_fore = inputs[Aircraft.Engine.WING_LOCATIONS][0][0]
 
         # TODO: high engine-count configuation.
@@ -349,7 +352,6 @@ def compute_partials(self, inputs, partials, discrete_inputs=None):
         # TODO temp using first engine, heterogeneous engines not supported
         num_eng = options.get_val(Aircraft.Engine.NUM_ENGINES)[0]
         num_wing_eng = options.get_val(Aircraft.Engine.NUM_WING_ENGINES)[0]
-
         y_eng_fore = inputs[Aircraft.Engine.WING_LOCATIONS][0][0]
         y_eng_aft = 0
 
diff --git a/aviary/subsystems/mass/flops_based/test/test_wing_common.py b/aviary/subsystems/mass/flops_based/test/test_wing_common.py
index 69134fa47..cfa95cb4a 100644
--- a/aviary/subsystems/mass/flops_based/test/test_wing_common.py
+++ b/aviary/subsystems/mass/flops_based/test/test_wing_common.py
@@ -159,24 +159,27 @@ def test_case(self, case_name):
         flops_validation_test(
             prob,
             case_name,
-            input_keys=[Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
-                        Aircraft.Wing.BENDING_FACTOR,
-                        Aircraft.Wing.BENDING_MASS_SCALER,
-                        Aircraft.Wing.COMPOSITE_FRACTION,
-                        Aircraft.Wing.ENG_POD_INERTIA_FACTOR,
-                        Mission.Design.GROSS_MASS,
-                        Aircraft.Wing.LOAD_FRACTION,
-                        Aircraft.Wing.MISC_MASS,
-                        Aircraft.Wing.MISC_MASS_SCALER,
-                        Aircraft.Wing.SHEAR_CONTROL_MASS,
-                        Aircraft.Wing.SHEAR_CONTROL_MASS_SCALER,
-                        Aircraft.Wing.SPAN,
-                        Aircraft.Wing.SWEEP,
-                        Aircraft.Wing.ULTIMATE_LOAD_FACTOR,
-                        Aircraft.Wing.VAR_SWEEP_MASS_PENALTY],
-            output_keys=Aircraft.Wing.BENDING_MASS,
+            input_keys=[
+                Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+                Aircraft.Wing.BENDING_MATERIAL_FACTOR,
+                Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER,
+                Aircraft.Wing.COMPOSITE_FRACTION,
+                Aircraft.Wing.ENG_POD_INERTIA_FACTOR,
+                Mission.Design.GROSS_MASS,
+                Aircraft.Wing.LOAD_FRACTION,
+                Aircraft.Wing.MISC_MASS,
+                Aircraft.Wing.MISC_MASS_SCALER,
+                Aircraft.Wing.SHEAR_CONTROL_MASS,
+                Aircraft.Wing.SHEAR_CONTROL_MASS_SCALER,
+                Aircraft.Wing.SPAN,
+                Aircraft.Wing.SWEEP,
+                Aircraft.Wing.ULTIMATE_LOAD_FACTOR,
+                Aircraft.Wing.VAR_SWEEP_MASS_PENALTY,
+            ],
+            output_keys=Aircraft.Wing.BENDING_MATERIAL_MASS,
             atol=1e-11,
-            rtol=1e-11)
+            rtol=1e-11,
+        )
 
     def test_IO(self):
         assert_match_varnames(self.prob.model)
@@ -206,7 +209,7 @@ def test_case(self):
         )
         prob.setup(check=False, force_alloc_complex=True)
         prob.set_val(Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR, 0.333, 'unitless')
-        prob.set_val(Aircraft.Wing.BENDING_FACTOR, 10, 'unitless')
+        prob.set_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR, 10, 'unitless')
         prob.set_val(Aircraft.Wing.COMPOSITE_FRACTION, 0.333, 'unitless')
         prob.set_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR, 1, 'unitless')
         prob.set_val(Mission.Design.GROSS_MASS, 100000, 'lbm')
diff --git a/aviary/subsystems/mass/flops_based/test/test_wing_detailed.py b/aviary/subsystems/mass/flops_based/test/test_wing_detailed.py
index fbb486d83..1b10c1e39 100644
--- a/aviary/subsystems/mass/flops_based/test/test_wing_detailed.py
+++ b/aviary/subsystems/mass/flops_based/test/test_wing_detailed.py
@@ -25,8 +25,12 @@ def setUp(self):
         self.prob = om.Problem()
 
     # Skip model that doesn't use detailed wing.
-    @parameterized.expand(get_flops_case_names(omit=['LargeSingleAisle2FLOPS', 'LargeSingleAisle2FLOPSalt']),
-                          name_func=print_case)
+    @parameterized.expand(
+        get_flops_case_names(
+            omit=[
+                'LargeSingleAisle2FLOPS',
+                'LargeSingleAisle2FLOPSalt']),
+        name_func=print_case)
     def test_case(self, case_name):
 
         prob = self.prob
@@ -43,23 +47,28 @@ def test_case(self, case_name):
         flops_validation_test(
             prob,
             case_name,
-            input_keys=[Aircraft.Wing.LOAD_PATH_SWEEP_DIST,
-                        Aircraft.Wing.THICKNESS_TO_CHORD_DIST,
-                        Aircraft.Wing.CHORD_PER_SEMISPAN_DIST,
-                        Mission.Design.GROSS_MASS,
-                        Aircraft.Engine.POD_MASS,
-                        Aircraft.Wing.ASPECT_RATIO,
-                        Aircraft.Wing.ASPECT_RATIO_REF,
-                        Aircraft.Wing.STRUT_BRACING_FACTOR,
-                        Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
-                        Aircraft.Engine.WING_LOCATIONS,
-                        Aircraft.Wing.THICKNESS_TO_CHORD,
-                        Aircraft.Wing.THICKNESS_TO_CHORD_REF],
-            output_keys=[Aircraft.Wing.BENDING_FACTOR,
-                         Aircraft.Wing.ENG_POD_INERTIA_FACTOR],
+            input_keys=[
+                Aircraft.Wing.LOAD_PATH_SWEEP_DIST,
+                Aircraft.Wing.THICKNESS_TO_CHORD_DIST,
+                Aircraft.Wing.CHORD_PER_SEMISPAN_DIST,
+                Mission.Design.GROSS_MASS,
+                Aircraft.Engine.POD_MASS,
+                Aircraft.Wing.ASPECT_RATIO,
+                Aircraft.Wing.ASPECT_RATIO_REF,
+                Aircraft.Wing.STRUT_BRACING_FACTOR,
+                Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+                Aircraft.Engine.WING_LOCATIONS,
+                Aircraft.Wing.THICKNESS_TO_CHORD,
+                Aircraft.Wing.THICKNESS_TO_CHORD_REF,
+            ],
+            output_keys=[
+                Aircraft.Wing.BENDING_MATERIAL_FACTOR,
+                Aircraft.Wing.ENG_POD_INERTIA_FACTOR,
+            ],
             method='fd',
             atol=1e-3,
-            rtol=1e-5)
+            rtol=1e-5,
+        )
 
     def test_case_multiengine(self):
         prob = self.prob
@@ -111,20 +120,141 @@ def test_case_multiengine(self):
 
         prob.run_model()
 
-        bending_factor = prob.get_val(Aircraft.Wing.BENDING_FACTOR)
+        BENDING_MATERIAL_FACTOR = prob.get_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR)
         pod_inertia = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
 
-        # manual computation of expected thrust reverser mass
-        bending_factor_expected = 11.59165669761
+        BENDING_MATERIAL_FACTOR_expected = 11.59165669761
         # 0.9600334354133278 if the factors are additive
         pod_inertia_expected = 0.9604608395586276
-        assert_near_equal(bending_factor, bending_factor_expected, tolerance=1e-10)
+        assert_near_equal(
+            BENDING_MATERIAL_FACTOR, BENDING_MATERIAL_FACTOR_expected, tolerance=1e-10
+        )
         assert_near_equal(pod_inertia, pod_inertia_expected, tolerance=1e-10)
 
         partial_data = prob.check_partials(
-            out_stream=None, compact_print=True, show_only_incorrect=True, form='central', method="fd")
+            out_stream=None,
+            compact_print=True,
+            show_only_incorrect=True,
+            form='central',
+            method="fd")
         assert_check_partials(partial_data, atol=1e-5, rtol=1e-5)
 
+    def test_case_fuselage_engines(self):
+        prob = self.prob
+
+        aviary_options = get_flops_inputs('LargeSingleAisle1FLOPS')
+
+        engine_options = AviaryValues()
+        engine_options.set_val(Settings.VERBOSITY, 0)
+        engine_options.set_val(Aircraft.Engine.DATA_FILE,
+                               'models/engines/turbofan_28k.deck')
+        engine_options.set_val(Aircraft.Engine.NUM_ENGINES, 2)
+        engine_options.set_val(Aircraft.Engine.NUM_WING_ENGINES, 0)
+        engine_options.set_val(Aircraft.Engine.NUM_FUSELAGE_ENGINES, 2)
+        engineModel = EngineDeck(options=engine_options)
+
+        preprocess_propulsion(aviary_options, [engineModel])
+
+        prob.model.add_subsystem('detailed_wing', DetailedWingBendingFact(
+            aviary_options=aviary_options), promotes=['*'])
+
+        prob.setup(force_alloc_complex=True)
+
+        input_keys = [Aircraft.Wing.LOAD_PATH_SWEEP_DIST,
+                      Aircraft.Wing.THICKNESS_TO_CHORD_DIST,
+                      Aircraft.Wing.CHORD_PER_SEMISPAN_DIST,
+                      Mission.Design.GROSS_MASS,
+                      Aircraft.Wing.ASPECT_RATIO,
+                      Aircraft.Wing.ASPECT_RATIO_REF,
+                      Aircraft.Wing.STRUT_BRACING_FACTOR,
+                      Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+                      Aircraft.Wing.THICKNESS_TO_CHORD,
+                      Aircraft.Wing.THICKNESS_TO_CHORD_REF]
+
+        for key in input_keys:
+            val, units = aviary_options.get_item(key)
+            prob.set_val(key, val, units)
+
+        prob.set_val(Aircraft.Engine.POD_MASS, np.array([0]), units='lbm')
+
+        wing_location = np.zeros(0)
+        wing_location = np.append(wing_location, [0.0])
+
+        prob.set_val(Aircraft.Engine.WING_LOCATIONS, wing_location)
+
+        prob.run_model()
+
+        BENDING_MATERIAL_FACTOR = prob.get_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR)
+        pod_inertia = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
+
+        BENDING_MATERIAL_FACTOR_expected = 11.59165669761
+        pod_inertia_expected = 0.84
+        assert_near_equal(
+            BENDING_MATERIAL_FACTOR, BENDING_MATERIAL_FACTOR_expected, tolerance=1e-10
+        )
+        assert_near_equal(pod_inertia, pod_inertia_expected, tolerance=1e-10)
+
+    def test_case_fuselage_multiengine(self):
+        prob = self.prob
+
+        aviary_options = get_flops_inputs('LargeSingleAisle1FLOPS')
+
+        engine_options = AviaryValues()
+        engine_options.set_val(Aircraft.Engine.DATA_FILE,
+                               'models/engines/turbofan_28k.deck')
+        engine_options.set_val(Aircraft.Engine.NUM_ENGINES, 2)
+        engine_options.set_val(Aircraft.Engine.NUM_WING_ENGINES, 2)
+        engine_options.set_val(Aircraft.Engine.NUM_FUSELAGE_ENGINES, 0)
+        engineModel1 = EngineDeck(options=engine_options)
+
+        engine_options.set_val(Aircraft.Engine.DATA_FILE,
+                               'models/engines/turbofan_22k.deck')
+        engine_options.set_val(Aircraft.Engine.NUM_ENGINES, 2)
+        engine_options.set_val(Aircraft.Engine.NUM_WING_ENGINES, 0)
+        engine_options.set_val(Aircraft.Engine.NUM_FUSELAGE_ENGINES, 2)
+        engineModel2 = EngineDeck(options=engine_options)
+
+        preprocess_propulsion(aviary_options, [engineModel1, engineModel2])
+
+        prob.model.add_subsystem('detailed_wing', DetailedWingBendingFact(
+            aviary_options=aviary_options), promotes=['*'])
+
+        prob.setup(force_alloc_complex=True)
+
+        input_keys = [Aircraft.Wing.LOAD_PATH_SWEEP_DIST,
+                      Aircraft.Wing.THICKNESS_TO_CHORD_DIST,
+                      Aircraft.Wing.CHORD_PER_SEMISPAN_DIST,
+                      Mission.Design.GROSS_MASS,
+                      Aircraft.Wing.ASPECT_RATIO,
+                      Aircraft.Wing.ASPECT_RATIO_REF,
+                      Aircraft.Wing.STRUT_BRACING_FACTOR,
+                      Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+                      Aircraft.Wing.THICKNESS_TO_CHORD,
+                      Aircraft.Wing.THICKNESS_TO_CHORD_REF]
+
+        for key in input_keys:
+            val, units = aviary_options.get_item(key)
+            prob.set_val(key, val, units)
+
+        prob.set_val(Aircraft.Engine.POD_MASS, np.array([1130, 0]), units='lbm')
+
+        wing_locations = np.zeros(0)
+        wing_locations = np.append(wing_locations, [0.5])
+
+        prob.set_val(Aircraft.Engine.WING_LOCATIONS, wing_locations)
+
+        prob.run_model()
+
+        BENDING_MATERIAL_FACTOR = prob.get_val(Aircraft.Wing.BENDING_MATERIAL_FACTOR)
+        pod_inertia = prob.get_val(Aircraft.Wing.ENG_POD_INERTIA_FACTOR)
+
+        BENDING_MATERIAL_FACTOR_expected = 11.59165669761
+        pod_inertia_expected = 0.84
+        assert_near_equal(
+            BENDING_MATERIAL_FACTOR, BENDING_MATERIAL_FACTOR_expected, tolerance=1e-10
+        )
+        assert_near_equal(pod_inertia, pod_inertia_expected, tolerance=1e-10)
+
     def test_extreme_engine_loc(self):
 
         aviary_options = get_flops_inputs('LargeSingleAisle1FLOPS')
@@ -190,4 +320,4 @@ def test_IO(self):
     # test = DetailedWingBendingTest()
     # test.setUp()
     # test.test_case(case_name='LargeSingleAisle1FLOPS')
-    # test.test_extreme_engine_loc()
+    # test.test_case_multiengine()
diff --git a/aviary/subsystems/mass/flops_based/test/test_wing_simple.py b/aviary/subsystems/mass/flops_based/test/test_wing_simple.py
index da5d4148b..4485c5954 100644
--- a/aviary/subsystems/mass/flops_based/test/test_wing_simple.py
+++ b/aviary/subsystems/mass/flops_based/test/test_wing_simple.py
@@ -36,18 +36,23 @@ def test_case(self, case_name):
         flops_validation_test(
             prob,
             case_name,
-            input_keys=[Aircraft.Wing.AREA,
-                        Aircraft.Wing.SPAN,
-                        Aircraft.Wing.TAPER_RATIO,
-                        Aircraft.Wing.THICKNESS_TO_CHORD,
-                        Aircraft.Wing.STRUT_BRACING_FACTOR,
-                        Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
-                        Aircraft.Wing.ASPECT_RATIO,
-                        Aircraft.Wing.SWEEP],
-            output_keys=[Aircraft.Wing.BENDING_FACTOR,
-                         Aircraft.Wing.ENG_POD_INERTIA_FACTOR],
+            input_keys=[
+                Aircraft.Wing.AREA,
+                Aircraft.Wing.SPAN,
+                Aircraft.Wing.TAPER_RATIO,
+                Aircraft.Wing.THICKNESS_TO_CHORD,
+                Aircraft.Wing.STRUT_BRACING_FACTOR,
+                Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+                Aircraft.Wing.ASPECT_RATIO,
+                Aircraft.Wing.SWEEP,
+            ],
+            output_keys=[
+                Aircraft.Wing.BENDING_MATERIAL_FACTOR,
+                Aircraft.Wing.ENG_POD_INERTIA_FACTOR,
+            ],
             atol=1e-11,
-            rtol=1e-11)
+            rtol=1e-11,
+        )
 
     def test_IO(self):
         assert_match_varnames(self.prob.model)
diff --git a/aviary/subsystems/mass/flops_based/wing_common.py b/aviary/subsystems/mass/flops_based/wing_common.py
index e810dad45..f673e29f3 100644
--- a/aviary/subsystems/mass/flops_based/wing_common.py
+++ b/aviary/subsystems/mass/flops_based/wing_common.py
@@ -21,8 +21,8 @@ def initialize(self):
     def setup(self):
         add_aviary_input(self, Mission.Design.GROSS_MASS, val=0.0)
         add_aviary_input(self, Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR, val=0.0)
-        add_aviary_input(self, Aircraft.Wing.BENDING_FACTOR, val=0.0)
-        add_aviary_input(self, Aircraft.Wing.BENDING_MASS_SCALER, val=1.0)
+        add_aviary_input(self, Aircraft.Wing.BENDING_MATERIAL_FACTOR, val=0.0)
+        add_aviary_input(self, Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER, val=1.0)
         add_aviary_input(self, Aircraft.Wing.COMPOSITE_FRACTION, val=0.0)
         add_aviary_input(self, Aircraft.Wing.ENG_POD_INERTIA_FACTOR, val=0.0)
         add_aviary_input(self, Aircraft.Wing.LOAD_FRACTION, val=0.0)
@@ -35,7 +35,7 @@ def setup(self):
         add_aviary_input(self, Aircraft.Wing.ULTIMATE_LOAD_FACTOR, val=3.75)
         add_aviary_input(self, Aircraft.Wing.VAR_SWEEP_MASS_PENALTY, val=0.0)
 
-        add_aviary_output(self, Aircraft.Wing.BENDING_MASS, val=0.0)
+        add_aviary_output(self, Aircraft.Wing.BENDING_MATERIAL_MASS, val=0.0)
 
         self.A1 = 8.80
         self.A2 = 6.25
@@ -44,7 +44,7 @@ def setup_partials(self):
         self.declare_partials("*", "*")
 
     def compute(self, inputs, outputs):
-        bt = inputs[Aircraft.Wing.BENDING_FACTOR]
+        bt = inputs[Aircraft.Wing.BENDING_MATERIAL_FACTOR]
         ulf = inputs[Aircraft.Wing.ULTIMATE_LOAD_FACTOR]
         span = inputs[Aircraft.Wing.SPAN]
         comp_frac = inputs[Aircraft.Wing.COMPOSITE_FRACTION]
@@ -54,7 +54,7 @@ def compute(self, inputs, outputs):
         sweep = inputs[Aircraft.Wing.SWEEP]
         gross_weight = inputs[Mission.Design.GROSS_MASS] * GRAV_ENGLISH_LBM
         CAYE = inputs[Aircraft.Wing.ENG_POD_INERTIA_FACTOR]
-        scaler = inputs[Aircraft.Wing.BENDING_MASS_SCALER]
+        scaler = inputs[Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER]
 
         num_fuse = self.options['aviary_options'].get_val(
             Aircraft.Fuselage.NUM_FUSELAGES)
@@ -72,12 +72,14 @@ def compute(self, inputs, outputs):
         W1NIR = self.A1 * bt * (1.0 + (self.A2 / span)**0.5) * ulf * span * \
             (1.0 - 0.4 * comp_frac) * (1.0 - 0.1 * faert) * cayf * vfact * pctl * 1.0e-6
 
-        outputs[Aircraft.Wing.BENDING_MASS] = (
-            (gross_weight * CAYE * W1NIR + W2 + W3) / (1.0 + W1NIR) - W2 - W3) \
-            * scaler / GRAV_ENGLISH_LBM
+        outputs[Aircraft.Wing.BENDING_MATERIAL_MASS] = (
+            ((gross_weight * CAYE * W1NIR + W2 + W3) / (1.0 + W1NIR) - W2 - W3)
+            * scaler
+            / GRAV_ENGLISH_LBM
+        )
 
     def compute_partials(self, inputs, J):
-        bt = inputs[Aircraft.Wing.BENDING_FACTOR]
+        bt = inputs[Aircraft.Wing.BENDING_MATERIAL_FACTOR]
         ulf = inputs[Aircraft.Wing.ULTIMATE_LOAD_FACTOR]
         span = inputs[Aircraft.Wing.SPAN]
         comp_frac = inputs[Aircraft.Wing.COMPOSITE_FRACTION]
@@ -91,7 +93,7 @@ def compute_partials(self, inputs, J):
         W3 = inputs[Aircraft.Wing.MISC_MASS] * GRAV_ENGLISH_LBM
         W2scale = inputs[Aircraft.Wing.SHEAR_CONTROL_MASS_SCALER]
         W3scale = inputs[Aircraft.Wing.MISC_MASS_SCALER]
-        scaler = inputs[Aircraft.Wing.BENDING_MASS_SCALER]
+        scaler = inputs[Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER]
 
         num_fuse = self.options['aviary_options'].get_val(
             Aircraft.Fuselage.NUM_FUSELAGES)
@@ -136,50 +138,69 @@ def compute_partials(self, inputs, J):
         fact2 = 1.0 / (1.0 + W1NIR)
         dbend_w1nir = scaler * (gross_weight * CAYE * fact2 - fact1 * fact2**2)
 
-        J[Aircraft.Wing.BENDING_MASS, Mission.Design.GROSS_MASS] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Mission.Design.GROSS_MASS] = (
             CAYE * W1NIR * fact2 * scaler
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.ENG_POD_INERTIA_FACTOR] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.ENG_POD_INERTIA_FACTOR] = (
             gross_weight * W1NIR * fact2 * scaler / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.SHEAR_CONTROL_MASS] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.SHEAR_CONTROL_MASS] = (
             (fact2 - 1.0) * scaler / W2scale
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.SHEAR_CONTROL_MASS_SCALER] = \
-            -(fact2 - 1.0) * scaler * W2 / W2scale ** 2 / GRAV_ENGLISH_LBM
+        J[
+            Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.SHEAR_CONTROL_MASS_SCALER
+        ] = (-(fact2 - 1.0) * scaler * W2 / W2scale**2 / GRAV_ENGLISH_LBM)
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.MISC_MASS] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.MISC_MASS] = (
             (fact2 - 1.0) * scaler / W3scale
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.MISC_MASS_SCALER] = \
-            -(fact2 - 1.0) * scaler * W3 / W3scale ** 2 / GRAV_ENGLISH_LBM
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.MISC_MASS_SCALER] = (
+            -(fact2 - 1.0) * scaler * W3 / W3scale**2 / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.BENDING_MASS_SCALER] = \
-            (fact1 * fact2 - W2/W2scale - W3/W3scale) / GRAV_ENGLISH_LBM
+        J[
+            Aircraft.Wing.BENDING_MATERIAL_MASS,
+            Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER,
+        ] = (fact1 * fact2 - W2 / W2scale - W3 / W3scale) / GRAV_ENGLISH_LBM
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.BENDING_FACTOR] = \
-            dbend_w1nir * dW1NIR_bt / GRAV_ENGLISH_LBM
+        J[
+            Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.BENDING_MATERIAL_FACTOR
+        ] = (dbend_w1nir * dW1NIR_bt / GRAV_ENGLISH_LBM)
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.ULTIMATE_LOAD_FACTOR] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.ULTIMATE_LOAD_FACTOR] = (
             dbend_w1nir * dW1NIR_ulf / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.LOAD_FRACTION] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.LOAD_FRACTION] = (
             dbend_w1nir * dW1NIR_pctl / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.COMPOSITE_FRACTION] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.COMPOSITE_FRACTION] = (
             dbend_w1nir * dW1NIR_compfrac / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR] = \
+        J[
+            Aircraft.Wing.BENDING_MATERIAL_MASS,
+            Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+        ] = (
             dbend_w1nir * dW1NIR_faert / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.VAR_SWEEP_MASS_PENALTY] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.VAR_SWEEP_MASS_PENALTY] = (
             dbend_w1nir * dW1NIR_varswp / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.SWEEP] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.SWEEP] = (
             dbend_w1nir * dW1NIR_sweep / GRAV_ENGLISH_LBM
+        )
 
-        J[Aircraft.Wing.BENDING_MASS, Aircraft.Wing.SPAN] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_MASS, Aircraft.Wing.SPAN] = (
             dbend_w1nir * dW1NIR_span / GRAV_ENGLISH_LBM
+        )
 
 
 class WingShearControlMass(om.ExplicitComponent):
@@ -336,7 +357,7 @@ def initialize(self):
             desc='collection of Aircraft/Mission specific options')
 
     def setup(self):
-        add_aviary_input(self, Aircraft.Wing.BENDING_MASS, val=0.0)
+        add_aviary_input(self, Aircraft.Wing.BENDING_MATERIAL_MASS, val=0.0)
 
         add_aviary_input(self, Aircraft.Wing.SHEAR_CONTROL_MASS, val=0.0)
 
@@ -352,7 +373,7 @@ def setup_partials(self):
         self.declare_partials("*", "*")
 
     def compute(self, inputs, outputs):
-        m1 = inputs[Aircraft.Wing.BENDING_MASS]
+        m1 = inputs[Aircraft.Wing.BENDING_MATERIAL_MASS]
         m2 = inputs[Aircraft.Wing.SHEAR_CONTROL_MASS]
         m3 = inputs[Aircraft.Wing.MISC_MASS]
         m4 = inputs[Aircraft.Wing.BWB_AFTBODY_MASS]
@@ -361,13 +382,13 @@ def compute(self, inputs, outputs):
         outputs[Aircraft.Wing.MASS] = (m1 + m2 + m3 + m4) * m_scaler
 
     def compute_partials(self, inputs, J):
-        m1 = inputs[Aircraft.Wing.BENDING_MASS]
+        m1 = inputs[Aircraft.Wing.BENDING_MATERIAL_MASS]
         m2 = inputs[Aircraft.Wing.SHEAR_CONTROL_MASS]
         m3 = inputs[Aircraft.Wing.MISC_MASS]
         m4 = inputs[Aircraft.Wing.BWB_AFTBODY_MASS]
         m_scaler = inputs[Aircraft.Wing.MASS_SCALER]
 
-        J[Aircraft.Wing.MASS, Aircraft.Wing.BENDING_MASS] = m_scaler
+        J[Aircraft.Wing.MASS, Aircraft.Wing.BENDING_MATERIAL_MASS] = m_scaler
         J[Aircraft.Wing.MASS, Aircraft.Wing.SHEAR_CONTROL_MASS] = m_scaler
         J[Aircraft.Wing.MASS, Aircraft.Wing.MISC_MASS] = m_scaler
         J[Aircraft.Wing.MASS, Aircraft.Wing.BWB_AFTBODY_MASS] = m_scaler
diff --git a/aviary/subsystems/mass/flops_based/wing_detailed.py b/aviary/subsystems/mass/flops_based/wing_detailed.py
index 7469ce55b..7763e50b9 100644
--- a/aviary/subsystems/mass/flops_based/wing_detailed.py
+++ b/aviary/subsystems/mass/flops_based/wing_detailed.py
@@ -16,15 +16,18 @@ class DetailedWingBendingFact(om.ExplicitComponent):
 
     def initialize(self):
         self.options.declare(
-            'aviary_options', types=AviaryValues,
-            desc='collection of Aircraft/Mission specific options')
+            'aviary_options',
+            types=AviaryValues,
+            desc='collection of Aircraft/Mission specific options',
+        )
 
     def setup(self):
         aviary_options: AviaryValues = self.options['aviary_options']
         input_station_dist = aviary_options.get_val(Aircraft.Wing.INPUT_STATION_DIST)
         num_input_stations = len(input_station_dist)
         total_num_wing_engines = aviary_options.get_val(
-            Aircraft.Propulsion.TOTAL_NUM_WING_ENGINES)
+            Aircraft.Propulsion.TOTAL_NUM_WING_ENGINES
+        )
         num_engine_type = len(aviary_options.get_val(Aircraft.Engine.NUM_ENGINES))
 
         # wing locations are different for each engine type - ragged array!
@@ -33,14 +36,23 @@ def setup(self):
         # wing_location_default = np.empty(num_engine_type, object)
         # wing_location_default[:] = [np.array([0]*int(num)) for num in num_wing_engines/2]
 
-        add_aviary_input(self, Aircraft.Wing.LOAD_PATH_SWEEP_DIST,
-                         val=np.zeros(num_input_stations - 1))
+        add_aviary_input(
+            self,
+            Aircraft.Wing.LOAD_PATH_SWEEP_DIST,
+            val=np.zeros(num_input_stations - 1),
+        )
 
-        add_aviary_input(self, Aircraft.Wing.THICKNESS_TO_CHORD_DIST,
-                         val=np.zeros(num_input_stations))
+        add_aviary_input(
+            self,
+            Aircraft.Wing.THICKNESS_TO_CHORD_DIST,
+            val=np.zeros(num_input_stations),
+        )
 
-        add_aviary_input(self, Aircraft.Wing.CHORD_PER_SEMISPAN_DIST,
-                         val=np.zeros(num_input_stations))
+        add_aviary_input(
+            self,
+            Aircraft.Wing.CHORD_PER_SEMISPAN_DIST,
+            val=np.zeros(num_input_stations),
+        )
 
         add_aviary_input(self, Mission.Design.GROSS_MASS, val=0.0)
 
@@ -54,14 +66,20 @@ def setup(self):
 
         add_aviary_input(self, Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR, val=0.0)
 
-        add_aviary_input(self, Aircraft.Engine.WING_LOCATIONS,
-                         val=np.zeros(int(total_num_wing_engines/2)))
+        if total_num_wing_engines > 0:
+            add_aviary_input(
+                self,
+                Aircraft.Engine.WING_LOCATIONS,
+                val=np.zeros(int(total_num_wing_engines / 2)),
+            )
+        else:
+            add_aviary_input(self, Aircraft.Engine.WING_LOCATIONS, val=[[0.0]])
 
         add_aviary_input(self, Aircraft.Wing.THICKNESS_TO_CHORD, val=0.0)
 
         add_aviary_input(self, Aircraft.Wing.THICKNESS_TO_CHORD_REF, val=0.0)
 
-        add_aviary_output(self, Aircraft.Wing.BENDING_FACTOR, val=0.0)
+        add_aviary_output(self, Aircraft.Wing.BENDING_MATERIAL_FACTOR, val=0.0)
 
         add_aviary_output(self, Aircraft.Wing.ENG_POD_INERTIA_FACTOR, val=0.0)
 
@@ -73,8 +91,9 @@ def compute(self, inputs, outputs):
         aviary_options: AviaryValues = self.options['aviary_options']
         input_station_dist = aviary_options.get_val(Aircraft.Wing.INPUT_STATION_DIST)
         inp_stations = np.array(input_station_dist)
-        num_integration_stations = \
-            aviary_options.get_val(Aircraft.Wing.NUM_INTEGRATION_STATIONS)
+        num_integration_stations = aviary_options.get_val(
+            Aircraft.Wing.NUM_INTEGRATION_STATIONS
+        )
         num_wing_engines = aviary_options.get_val(Aircraft.Engine.NUM_WING_ENGINES)
         num_engine_type = len(num_wing_engines)
 
@@ -85,8 +104,9 @@ def compute(self, inputs, outputs):
         # 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 = \
-            aviary_options.get_val(Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL)
+        load_distribution_factor = aviary_options.get_val(
+            Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL
+        )
 
         load_path_sweep = inputs[Aircraft.Wing.LOAD_PATH_SWEEP_DIST]
         thickness_to_chord = inputs[Aircraft.Wing.THICKNESS_TO_CHORD_DIST]
@@ -121,53 +141,77 @@ def compute(self, inputs, outputs):
             per_section = int(stations_per_section[i])
             integration_stations = np.append(
                 integration_stations,
-                np.linspace(inp_stations[i], val, per_section, endpoint=endpoint))
-            sweep_int_stations = np.append(sweep_int_stations,
-                                           load_path_sweep[i] * np.ones(per_section))
+                np.linspace(inp_stations[i], val, per_section, endpoint=endpoint),
+            )
+            sweep_int_stations = np.append(
+                sweep_int_stations, load_path_sweep[i] * np.ones(per_section)
+            )
 
         dy = np.diff(integration_stations)
-        avg_sweep = np.sum((dy[1:] + 2.0 * integration_stations[1:-1]) * dy[1:] *
-                           sweep_int_stations[1:-1])
+        avg_sweep = np.sum(
+            (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)
+            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 {Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL}, it must be "1", "2", or "3".')
+                f'{load_distribution_factor} is not a valid value for '
+                f'{Aircraft.Wing.LOAD_DISTRIBUTION_CONTROL}, it must be "1", "2", or "3".'
+            )
 
-        chord_interp = InterpND(method='slinear', points=(inp_stations),
-                                x_interp=integration_stations)
+        chord_interp = InterpND(
+            method='slinear', points=(inp_stations), x_interp=integration_stations
+        )
         chord_int_stations = chord_interp.evaluate_spline(
-            chord, compute_derivative=False)
+            chord, compute_derivative=False
+        )
         if arref > 0.0:
             # Scale
             chord_int_stations *= arref / ar
 
-        del_load = dy * (
-            chord_int_stations[:-1] * (2*load_intensity[:-1] + load_intensity[1:]) +
-            chord_int_stations[1:] * (2*load_intensity[1:] + load_intensity[:-1])) / 6
+        del_load = (
+            dy
+            * (
+                chord_int_stations[:-1] * (2 * load_intensity[:-1] + load_intensity[1:])
+                + chord_int_stations[1:]
+                * (2 * load_intensity[1:] + load_intensity[:-1])
+            )
+            / 6
+        )
 
         el = np.sum(del_load)
 
-        del_moment = dy**2 * (
-            chord_int_stations[:-1] * (load_intensity[:-1]+load_intensity[1:]) +
-            chord_int_stations[1:] * (3*load_intensity[1:]+load_intensity[:-1])) / 12
+        del_moment = (
+            dy**2
+            * (
+                chord_int_stations[:-1] * (load_intensity[:-1] + load_intensity[1:])
+                + chord_int_stations[1:]
+                * (3 * load_intensity[1:] + load_intensity[:-1])
+            )
+            / 12
+        )
 
         load_path_length = np.flip(
-            np.append(np.zeros(1, chord.dtype), np.cumsum(np.flip(del_load)[:-1])))
-        csw = 1. / np.cos(sweep_int_stations[:-1] * np.pi/180.)
+            np.append(np.zeros(1, chord.dtype), np.cumsum(np.flip(del_load)[:-1]))
+        )
+        csw = 1.0 / np.cos(sweep_int_stations[:-1] * np.pi / 180.0)
         emi = (del_moment + dy * load_path_length) * csw
         # em = np.sum(emi)
 
-        tc_interp = InterpND(method='slinear', points=(inp_stations),
-                             x_interp=integration_stations)
+        tc_interp = InterpND(
+            method='slinear', points=(inp_stations), x_interp=integration_stations
+        )
         tc_int_stations = tc_interp.evaluate_spline(
-            thickness_to_chord, compute_derivative=False)
+            thickness_to_chord, compute_derivative=False
+        )
         if tcref > 0.0:
             tc_int_stations *= tc / tcref
 
@@ -181,16 +225,22 @@ def compute(self, inputs, outputs):
 
         btb = 4 * pm / el
 
-        sa = np.sin(avg_sweep * np.pi / 180.)
+        sa = np.sin(avg_sweep * np.pi / 180.0)
         if ar <= 5.0:
             caya = 0.0
         else:
             caya = ar - 5.0
 
-        bt = btb / (ar**(0.25*fstrt) * (1.0 + (0.5*faert - 0.16*fstrt)
-                    * sa**2 + 0.03*caya * (1.0-0.5*faert)*sa))
+        bt = btb / (
+            ar ** (0.25 * fstrt)
+            * (
+                1.0
+                + (0.5 * faert - 0.16 * fstrt) * sa**2
+                + 0.03 * caya * (1.0 - 0.5 * faert) * sa
+            )
+        )
 
-        outputs[Aircraft.Wing.BENDING_FACTOR] = bt
+        outputs[Aircraft.Wing.BENDING_MATERIAL_FACTOR] = bt
 
         inertia_factor = np.zeros(num_engine_type, dtype=chord.dtype)
         eel = np.zeros(len(dy) + 1, dtype=chord.dtype)
@@ -200,25 +250,29 @@ def compute(self, inputs, outputs):
         # i is the counter for which engine model we are checking
         for i in range(num_engine_type):
             # idx2 is the last index for the range of engines of this type
-            idx2 = idx + int(num_wing_engines[i]/2)
+            idx2 = idx + int(num_wing_engines[i] / 2)
+            if num_wing_engines[i] > 0:
+                # engine locations must be in order from wing root to tip
+                eng_loc = np.sort(engine_locations[idx:idx2])
 
-            eng_loc = engine_locations[idx:idx2][0]
+            else:
+                continue
 
-            if eng_loc <= integration_stations[0]:
+            if eng_loc[0] <= integration_stations[0]:
                 inertia_factor[i] = 1.0
 
-            elif eng_loc >= integration_stations[-1]:
+            elif eng_loc[0] >= integration_stations[-1]:
                 inertia_factor[i] = 0.84
 
             else:
                 eel[:] = 0.0
-                loc = np.where(integration_stations < eng_loc)[0]
+                # Find all points on integration station before first engine
+                loc = np.where(integration_stations < eng_loc[0])[0]
                 eel[loc] = 1.0
 
                 delme = dy * eel[1:]
 
-                delme[loc[-1]] = engine_locations[idx:idx2][0] - \
-                    integration_stations[loc[-1]]
+                delme[loc[-1]] = eng_loc[0] - integration_stations[loc[-1]]
 
                 eem = delme * csw
                 eem = np.cumsum(eem[::-1])[::-1]
diff --git a/aviary/subsystems/mass/flops_based/wing_group.py b/aviary/subsystems/mass/flops_based/wing_group.py
index 9be1571f3..284934f8f 100644
--- a/aviary/subsystems/mass/flops_based/wing_group.py
+++ b/aviary/subsystems/mass/flops_based/wing_group.py
@@ -28,14 +28,20 @@ def setup(self):
                            promotes_inputs=['*'], promotes_outputs=['*'])
 
         if Aircraft.Wing.INPUT_STATION_DIST in aviary_options:
-            self.add_subsystem('wing_bending_factor',
-                               DetailedWingBendingFact(aviary_options=aviary_options),
-                               promotes_inputs=['*'], promotes_outputs=['*'])
+            self.add_subsystem(
+                'wing_BENDING_MATERIAL_FACTOR',
+                DetailedWingBendingFact(aviary_options=aviary_options),
+                promotes_inputs=['*'],
+                promotes_outputs=['*'],
+            )
 
         else:
-            self.add_subsystem('wing_bending_factor',
-                               SimpleWingBendingFact(aviary_options=aviary_options),
-                               promotes_inputs=['*'], promotes_outputs=['*'])
+            self.add_subsystem(
+                'wing_BENDING_MATERIAL_FACTOR',
+                SimpleWingBendingFact(aviary_options=aviary_options),
+                promotes_inputs=['*'],
+                promotes_outputs=['*'],
+            )
 
         self.add_subsystem('wing_misc',
                            WingMiscMass(aviary_options=aviary_options),
diff --git a/aviary/subsystems/mass/flops_based/wing_simple.py b/aviary/subsystems/mass/flops_based/wing_simple.py
index 180c8f3de..d41f5d321 100644
--- a/aviary/subsystems/mass/flops_based/wing_simple.py
+++ b/aviary/subsystems/mass/flops_based/wing_simple.py
@@ -34,20 +34,24 @@ def setup(self):
 
         add_aviary_input(self, Aircraft.Wing.SWEEP, val=0.0)
 
-        add_aviary_output(self, Aircraft.Wing.BENDING_FACTOR, val=0.0)
+        add_aviary_output(self, Aircraft.Wing.BENDING_MATERIAL_FACTOR, val=0.0)
 
         add_aviary_output(self, Aircraft.Wing.ENG_POD_INERTIA_FACTOR, val=0.0)
 
     def setup_partials(self):
-        self.declare_partials(of=Aircraft.Wing.BENDING_FACTOR,
-                              wrt=[Aircraft.Wing.STRUT_BRACING_FACTOR,
-                                   Aircraft.Wing.SPAN,
-                                   Aircraft.Wing.TAPER_RATIO,
-                                   Aircraft.Wing.AREA,
-                                   Aircraft.Wing.THICKNESS_TO_CHORD,
-                                   Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
-                                   Aircraft.Wing.ASPECT_RATIO,
-                                   Aircraft.Wing.SWEEP])
+        self.declare_partials(
+            of=Aircraft.Wing.BENDING_MATERIAL_FACTOR,
+            wrt=[
+                Aircraft.Wing.STRUT_BRACING_FACTOR,
+                Aircraft.Wing.SPAN,
+                Aircraft.Wing.TAPER_RATIO,
+                Aircraft.Wing.AREA,
+                Aircraft.Wing.THICKNESS_TO_CHORD,
+                Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+                Aircraft.Wing.ASPECT_RATIO,
+                Aircraft.Wing.SWEEP,
+            ],
+        )
 
     def compute(self, inputs, outputs):
         aviary_options: AviaryValues = self.options['aviary_options']
@@ -77,8 +81,9 @@ def compute(self, inputs, outputs):
 
         ems = 1.0 - 0.25 * fstrt
 
-        outputs[Aircraft.Wing.BENDING_FACTOR] = \
-            0.215 * (0.37 + 0.7 * tr) * (span**2 / area)**ems / (cayl * tca)
+        outputs[Aircraft.Wing.BENDING_MATERIAL_FACTOR] = (
+            0.215 * (0.37 + 0.7 * tr) * (span**2 / area) ** ems / (cayl * tca)
+        )
 
         outputs[Aircraft.Wing.ENG_POD_INERTIA_FACTOR] = 1.0 - 0.03 * num_wing_eng
 
@@ -134,26 +139,37 @@ def compute_partials(self, inputs, J):
         dbend_tr = 0.215 * 0.7 * term2 * term3
         dbend_cayl = -0.215 * term1 * term2 * tca * term3**2
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.STRUT_BRACING_FACTOR] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_FACTOR, Aircraft.Wing.STRUT_BRACING_FACTOR] = (
             dbend_exp + dbend_cayl * dcayl_fstrt
+        )
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.SPAN] = \
-            2.0 * 0.215 * term1 * ems * term2a**(ems - 1) * term3 * span / area
+        J[Aircraft.Wing.BENDING_MATERIAL_FACTOR, Aircraft.Wing.SPAN] = (
+            2.0 * 0.215 * term1 * ems * term2a ** (ems - 1) * term3 * span / area
+        )
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.TAPER_RATIO] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_FACTOR, Aircraft.Wing.TAPER_RATIO] = (
             dbend_tr + dbend_cayl * (dcayl_slam * dslam * dtlam_tr)
+        )
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.AREA] = \
-            -0.215 * term1 * ems * term2a**(ems - 1) * term3 * term2a / area
+        J[Aircraft.Wing.BENDING_MATERIAL_FACTOR, Aircraft.Wing.AREA] = (
+            -0.215 * term1 * ems * term2a ** (ems - 1) * term3 * term2a / area
+        )
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.THICKNESS_TO_CHORD] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_FACTOR, Aircraft.Wing.THICKNESS_TO_CHORD] = (
             -0.215 * term1 * term2 * cayl * term3**2
+        )
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR] = \
+        J[
+            Aircraft.Wing.BENDING_MATERIAL_FACTOR,
+            Aircraft.Wing.AEROELASTIC_TAILORING_FACTOR,
+        ] = (
             dbend_cayl * dcayl_faert
+        )
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.ASPECT_RATIO] = \
+        J[Aircraft.Wing.BENDING_MATERIAL_FACTOR, Aircraft.Wing.ASPECT_RATIO] = (
             dbend_cayl * (dcayl_ar + dcayl_slam * dslam * dtlam_ar)
+        )
 
-        J[Aircraft.Wing.BENDING_FACTOR, Aircraft.Wing.SWEEP] = \
-            dbend_cayl * (dcayl_slam * dslam * dtlam_sweep)
+        J[Aircraft.Wing.BENDING_MATERIAL_FACTOR, Aircraft.Wing.SWEEP] = dbend_cayl * (
+            dcayl_slam * dslam * dtlam_sweep
+        )
diff --git a/aviary/subsystems/mass/test/test_flops_mass_builder.py b/aviary/subsystems/mass/test/test_flops_mass_builder.py
index 2e9b0ee98..98b152507 100644
--- a/aviary/subsystems/mass/test/test_flops_mass_builder.py
+++ b/aviary/subsystems/mass/test/test_flops_mass_builder.py
@@ -26,9 +26,12 @@ def setUp(self):
             'test_core_mass', meta_data=BaseMetaData, code_origin=FLOPS)
         self.aviary_values = av.AviaryValues()
         self.aviary_values.set_val(Aircraft.Design.USE_ALT_MASS, False, units='unitless')
-        self.aviary_values.set_val(Aircraft.Engine.NUM_ENGINES, [1], units='unitless')
         self.aviary_values.set_val(
-            Aircraft.Engine.NUM_WING_ENGINES, [2], units='unitless')
+            Aircraft.Engine.NUM_ENGINES, np.array([1]), units='unitless'
+        )
+        self.aviary_values.set_val(
+            Aircraft.Engine.NUM_WING_ENGINES, np.array([2]), units='unitless'
+        )
 
 
 class TestFLOPSMassBuilderAltMass(av.TestSubsystemBuilderBase):
@@ -43,9 +46,12 @@ def setUp(self):
             'test_core_mass', meta_data=BaseMetaData, code_origin=FLOPS)
         self.aviary_values = av.AviaryValues()
         self.aviary_values.set_val(Aircraft.Design.USE_ALT_MASS, True, units='unitless')
-        self.aviary_values.set_val(Aircraft.Engine.NUM_ENGINES, [1], units='unitless')
         self.aviary_values.set_val(
-            Aircraft.Engine.NUM_WING_ENGINES, [2], units='unitless')
+            Aircraft.Engine.NUM_ENGINES, np.array([1]), units='unitless'
+        )
+        self.aviary_values.set_val(
+            Aircraft.Engine.NUM_WING_ENGINES, np.array([2]), units='unitless'
+        )
 
 
 if __name__ == '__main__':
diff --git a/aviary/variable_info/variable_meta_data.py b/aviary/variable_info/variable_meta_data.py
index 33b271aac..a87c69465 100644
--- a/aviary/variable_info/variable_meta_data.py
+++ b/aviary/variable_info/variable_meta_data.py
@@ -728,10 +728,11 @@
 add_meta_data(
     Aircraft.CrewPayload.Design.NUM_BUSINESS_CLASS,
     meta_data=_MetaData,
-    historical_name={"GASP": None,
-                     "FLOPS": 'WTIN.NPB',  # ['&DEFINE.WTIN.NPB', 'WTS.NPB'],
-                     "LEAPS1": 'aircraft.inputs.L0_crew_and_payload.business_class_count'
-                     },
+    historical_name={
+        "GASP": None,
+        "FLOPS": 'WTIN.NPB',  # ['&DEFINE.WTIN.NPB', 'WTS.NPB'],
+        "LEAPS1": 'aircraft.inputs.L0_crew_and_payload.business_class_count',
+    },
     units='unitless',
     desc='number of business class passengers that the aircraft is designed to accommodate',
     types=int,
@@ -742,10 +743,11 @@
 add_meta_data(
     Aircraft.CrewPayload.Design.NUM_FIRST_CLASS,
     meta_data=_MetaData,
-    historical_name={"GASP": None,
-                     "FLOPS": 'WTIN.NPF',  # ['&DEFINE.WTIN.NPF', 'WTS.NPF'],
-                     "LEAPS1": 'aircraft.inputs.L0_crew_and_payload.first_class_count'
-                     },
+    historical_name={
+        "GASP": None,
+        "FLOPS": 'WTIN.NPF',  # ['&DEFINE.WTIN.NPF', 'WTS.NPF'],
+        "LEAPS1": 'aircraft.inputs.L0_crew_and_payload.first_class_count',
+    },
     units='unitless',
     desc='number of first class passengers that the aircraft is designed to accommodate',
     types=int,
@@ -756,10 +758,11 @@
 add_meta_data(
     Aircraft.CrewPayload.Design.NUM_PASSENGERS,
     meta_data=_MetaData,
-    historical_name={"GASP": 'INGASP.PAX',  # number of passenger seats excluding crew
-                     "FLOPS": None,  # ['CSTDAT.NSV', '~WEIGHT.NPASS', '~WTSTAT.NPASS'],
-                     "LEAPS1": 'aircraft.outputs.L0_crew_and_payload.passenger_count'
-                     },
+    historical_name={
+        "GASP": 'INGASP.PAX',  # number of passenger seats excluding crew
+        "FLOPS": None,  # ['CSTDAT.NSV', '~WEIGHT.NPASS', '~WTSTAT.NPASS'],
+        "LEAPS1": 'aircraft.outputs.L0_crew_and_payload.passenger_count',
+    },
     units='unitless',
     desc='total number of passengers that the aircraft is designed to accommodate',
     option=True,
@@ -772,10 +775,11 @@
 add_meta_data(
     Aircraft.CrewPayload.Design.NUM_TOURIST_CLASS,
     meta_data=_MetaData,
-    historical_name={"GASP": None,
-                     "FLOPS": 'WTIN.NPT',  # ['&DEFINE.WTIN.NPT', 'WTS.NPT'],
-                     "LEAPS1": 'aircraft.inputs.L0_crew_and_payload.tourist_class_count'
-                     },
+    historical_name={
+        "GASP": None,
+        "FLOPS": 'WTIN.NPT',  # ['&DEFINE.WTIN.NPT', 'WTS.NPT'],
+        "LEAPS1": 'aircraft.inputs.L0_crew_and_payload.tourist_class_count',
+    },
     units='unitless',
     desc='number of tourist class passengers that the aircraft is designed to accommodate',
     types=int,
@@ -1023,10 +1027,7 @@
     meta_data=_MetaData,
     # note: this GASP variable does not include cargo, but it does include
     # passenger baggage
-    historical_name={"GASP": 'INGASP.WPL',
-                     "FLOPS": None,
-                     "LEAPS1": None
-                     },
+    historical_name={"GASP": 'INGASP.WPL', "FLOPS": None, "LEAPS1": None},
     units='lbm',
     desc='mass of passenger payload, including passengers, passenger baggage',
 )
@@ -1304,11 +1305,13 @@
 )
 
 add_meta_data(
-    Aircraft.Design.IJEFF, meta_data=_MetaData,
+    Aircraft.Design.IJEFF,
+    meta_data=_MetaData,
     historical_name={"GASP": 'INGASP.IJEFF', "FLOPS": None, "LEAPS1": None},
     desc="A flag used by Jeff V. Bowles to debug GASP code during his 53 years supporting the development of GASP. "
     "This flag is planted here to thank him for his hard work and dedication, Aviary wouldn't be what it is today "
-    "without his help.",)
+    "without his help.",
+)
 
 add_meta_data(
     Aircraft.Design.LAMINAR_FLOW_LOWER,
@@ -1439,15 +1442,18 @@
 )
 
 add_meta_data(
-    Aircraft.Design.RESERVE_FUEL_FRACTION, meta_data=_MetaData,
+    Aircraft.Design.RESERVE_FUEL_FRACTION,
+    meta_data=_MetaData,
     historical_name={"GASP": None, "FLOPS": None, "LEAPS1": None},
-    option=True, units="unitless",
+    option=True,
+    units="unitless",
     desc='required fuel reserves: given as a proportion of mission fuel. This value must be nonnegative. '
     'Mission fuel only includes normal phases and excludes reserve phases. '
     'If it is 0.5, the reserve fuel is half of the mission fuel (one third of the total fuel). Note '
     'it can be greater than 1. If it is 2, there would be twice as much reserve fuel as mission fuel '
     '(the total fuel carried would be 1/3 for the mission and 2/3 for the reserve)',
-    default_value=0,)
+    default_value=0,
+)
 
 add_meta_data(
     Aircraft.Design.SMOOTH_MASS_DISCONTINUITIES,
@@ -1935,7 +1941,8 @@
 # Global hybrid throttle is also False by default to account for parallel-hybrid engines
 # that can't operate at every power level at every condition due to other constraints
 add_meta_data(
-    Aircraft.Engine.GLOBAL_HYBRID_THROTTLE, meta_data=_MetaData,
+    Aircraft.Engine.GLOBAL_HYBRID_THROTTLE,
+    meta_data=_MetaData,
     historical_name={"GASP": None, "FLOPS": None, "LEAPS1": None},
     units='unitless',
     desc='Flag for engine decks if the range of provided hybrid throttles is consistent '
@@ -1946,7 +1953,7 @@
     'individually normalized from 0 to 1 independent of other points on the deck).',
     default_value=False,
     types=bool,
-    option=True
+    option=True,
 )
 
 # TODO Disabling global throttle ranges is preferred (therefore default) to prevent
@@ -1955,29 +1962,31 @@
 add_meta_data(
     Aircraft.Engine.GLOBAL_THROTTLE,
     meta_data=_MetaData,
-    historical_name={"GASP": None,
-                     "FLOPS": None,
-                     "LEAPS1": None
-                     },
+    historical_name={"GASP": None, "FLOPS": None, "LEAPS1": None},
     units='unitless',
     desc='Flag for engine decks if the range of provided throttles is consistent '
-         'across all flight conditions (e.g. the maximum throttle seen in the entire '
-         'deck is 1.0, but a given flight condition only goes to 0.9 -> GLOBAL_THROTTLE '
-         '= TRUE means the engine can be extrapolated out to 1.0 at that point. If '
-         "GLOBAL_THROTTLE is False, then each flight condition's throttle range is "
-         'individually normalized from 0 to 1 independent of other points on the deck).',
+    'across all flight conditions (e.g. the maximum throttle seen in the entire '
+    'deck is 1.0, but a given flight condition only goes to 0.9 -> GLOBAL_THROTTLE '
+    '= TRUE means the engine can be extrapolated out to 1.0 at that point. If '
+    "GLOBAL_THROTTLE is False, then each flight condition's throttle range is "
+    'individually normalized from 0 to 1 independent of other points on the deck).',
     default_value=False,
     types=bool,
-    option=True
+    option=True,
 )
 
 # TODO dependency on NTYE? Does this var need preprocessing? Can this mention be removed?
 add_meta_data(
-    Aircraft.Engine.HAS_PROPELLERS, meta_data=_MetaData,
+    Aircraft.Engine.HAS_PROPELLERS,
+    meta_data=_MetaData,
     historical_name={"GASP": None, "FLOPS": None, "LEAPS1": None},
-    option=True, units="unitless", default_value=False, types=bool,
+    option=True,
+    units="unitless",
+    default_value=False,
+    types=bool,
     desc='if True, the aircraft has propellers, otherwise aircraft is assumed to have no '
-    'propellers. In GASP this depended on NTYE',)
+    'propellers. In GASP this depended on NTYE',
+)
 
 add_meta_data(
     Aircraft.Engine.IGNORE_NEGATIVE_THRUST,
@@ -1996,11 +2005,16 @@
 )
 
 add_meta_data(
-    Aircraft.Engine.INTERPOLATION_METHOD, meta_data=_MetaData,
+    Aircraft.Engine.INTERPOLATION_METHOD,
+    meta_data=_MetaData,
     historical_name={"GASP": None, "FLOPS": None, "LEAPS1": None},
-    units="unitless", option=True, default_value='slinear', types=str,
+    units="unitless",
+    option=True,
+    default_value='slinear',
+    types=str,
     desc="method used for interpolation on an engine deck's data file, allowable values are "
-    'table methods from openmdao.components.interp_util.interp',)
+    'table methods from openmdao.components.interp_util.interp',
+)
 
 add_meta_data(
     Aircraft.Engine.MASS,
@@ -4987,7 +5001,7 @@
         "LEAPS1": 'aircraft.inputs.L0_detailed_wing.ref_aspect_ratio',
     },
     units='unitless',
-    desc='Reference aspect ratio, used for detailed wing bending.',
+    desc='Reference aspect ratio, used for detailed wing mass estimation.',
 )
 
 add_meta_data(
@@ -4999,7 +5013,7 @@
 )
 
 add_meta_data(
-    Aircraft.Wing.BENDING_FACTOR,
+    Aircraft.Wing.BENDING_MATERIAL_FACTOR,
     meta_data=_MetaData,
     historical_name={
         "GASP": None,
@@ -5007,13 +5021,14 @@
         "LEAPS1": 'aircraft.outputs.L0_wing.bending_material_factor',
     },
     units='unitless',
-    desc='wing bending factor',
+    desc='Wing bending material factor with sweep adjustment. Used to compute '
+    'Aircraft.Wing.BENDING_MATERIAL_MASS',
 )
 
 add_meta_data(
     # Note user override
-    #    - see also: Aircraft.Wing.BENDING_MASS_SCALER
-    Aircraft.Wing.BENDING_MASS,
+    #    - see also: Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER
+    Aircraft.Wing.BENDING_MATERIAL_MASS,
     meta_data=_MetaData,
     historical_name={
         "GASP": None,
@@ -5026,7 +5041,7 @@
 )
 
 add_meta_data(
-    Aircraft.Wing.BENDING_MASS_SCALER,
+    Aircraft.Wing.BENDING_MATERIAL_MASS_SCALER,
     meta_data=_MetaData,
     historical_name={
         "GASP": None,
@@ -5211,7 +5226,8 @@
         "LEAPS1": 'aircraft.outputs.L0_wing.engine_inertia_relief_factor',
     },
     units='unitless',
-    desc='engine inertia relief factor',
+    desc='Engine inertia relief factor for wingspan inboard of engine locations. Used '
+    'to compute Aircraft.Wing.BENDING_MATERIAL_MASS',
 )
 
 add_meta_data(
@@ -6021,7 +6037,7 @@
         "LEAPS1": 'aircraft.inputs.L0_detailed_wing.ref_thickness_to_chord_ratio',
     },
     units='unitless',
-    desc='Reference thickness-to-chord ratio, used for detailed wing bending.',
+    desc='Reference thickness-to-chord ratio, used for detailed wing mass estimation.',
     default_value=0.0,
 )
 
@@ -6161,12 +6177,9 @@
 add_meta_data(
     Dynamic.Atmosphere.KINEMATIC_VISCOSITY,
     meta_data=_MetaData,
-    historical_name={"GASP": 'XKV',
-                     "FLOPS": None,
-                     "LEAPS1": None
-                     },
+    historical_name={"GASP": 'XKV', "FLOPS": None, "LEAPS1": None},
     units='ft**2/s',
-    desc="Atmospheric kinematic viscosity at the vehicle's current flight condition"
+    desc="Atmospheric kinematic viscosity at the vehicle's current flight condition",
 )
 
 add_meta_data(
@@ -6245,23 +6258,17 @@
 add_meta_data(
     Dynamic.Mission.DISTANCE,
     meta_data=_MetaData,
-    historical_name={"GASP": None,
-                     "FLOPS": 'range',
-                     "LEAPS1": None
-                     },
+    historical_name={"GASP": None, "FLOPS": 'range', "LEAPS1": None},
     units='NM',
-    desc="The total distance the vehicle has traveled since brake release at the current time"
+    desc="The total distance the vehicle has traveled since brake release at the current time",
 )
 
 add_meta_data(
     Dynamic.Mission.DISTANCE_RATE,
     meta_data=_MetaData,
-    historical_name={"GASP": None,
-                     "FLOPS": 'range_rate',
-                     "LEAPS1": None
-                     },
+    historical_name={"GASP": None, "FLOPS": 'range_rate', "LEAPS1": None},
     units='NM/s',
-    desc="The rate at which the distance traveled is changing at the current time"
+    desc="The rate at which the distance traveled is changing at the current time",
 )
 
 add_meta_data(
@@ -7634,14 +7641,18 @@
 )
 
 add_meta_data(
-    Settings.VERBOSITY, meta_data=_MetaData,
+    Settings.VERBOSITY,
+    meta_data=_MetaData,
     historical_name={"GASP": None, "FLOPS": None, "LEAPS1": None},
     desc='Sets how much information Aviary outputs when run. Options include:'
     '0. QUIET: All output except errors are suppressed'
     '1. BRIEF: Only important information is output, in human-readable format'
     '2. VERBOSE: All user-relevant information is output, in human-readable format'
     '3. DEBUG: Any information can be outtputed, including warnings, intermediate calculations, etc., with no formatting requirement',
-    option=True, types=Verbosity, default_value=Verbosity.BRIEF,)
+    option=True,
+    types=Verbosity,
+    default_value=Verbosity.BRIEF,
+)
 
 # here we create a copy of the Aviary-core metadata. The reason for this
 # copy is that if we simply imported the Aviary _MetaData in all the
diff --git a/aviary/variable_info/variables.py b/aviary/variable_info/variables.py
index f9ad1c69e..9e8b2a709 100644
--- a/aviary/variable_info/variables.py
+++ b/aviary/variable_info/variables.py
@@ -508,11 +508,11 @@ class Wing:
         ASPECT_RATIO = 'aircraft:wing:aspect_ratio'
         ASPECT_RATIO_REF = 'aircraft:wing:aspect_ratio_reference'
         AVERAGE_CHORD = 'aircraft:wing:average_chord'
-        BENDING_FACTOR = 'aircraft:wing:bending_factor'
-        BENDING_MASS = 'aircraft:wing:bending_mass'
+        BENDING_MATERIAL_FACTOR = 'aircraft:wing:bending_material_factor'
+        BENDING_MATERIAL_MASS = 'aircraft:wing:bending_material_mass'
         # Not defined in metadata!
         # BENDING_MASS_NO_INERTIA = 'aircraft:wing:bending_mass_no_inertia'
-        BENDING_MASS_SCALER = 'aircraft:wing:bending_mass_scaler'
+        BENDING_MATERIAL_MASS_SCALER = 'aircraft:wing:bending_material_mass_scaler'
         BWB_AFTBODY_MASS = 'aircraft:wing:bwb_aft_body_mass'
         BWB_AFTBODY_MASS_SCALER = 'aircraft:wing:bwb_aft_body_mass_scaler'
         CENTER_CHORD = 'aircraft:wing:center_chord'