Update Hamilton Standard model based on Jeff's enhancement

aviary/subsystems/propulsion/hamilton_standard.py

@@ -88,6 +88,7 @@ def _biquad(T, i, xi, yi):
     T(i+3) = values of x in ascending order
     """
 
+    lmt = 0
     nx = int(T[i])
     ny = int(T[i+1])
     j1 = int(i + 2)
@@ -361,16 +362,16 @@ def _biquad(T, i, xi, yi):
 CL_arr = np.array([0.3, 0.4, 0.5, 0.6, 0.7, 0.8])
 CP_CLi_table = np.array([
     [0.0114, 0.0294, .0491, .0698, .0913, .1486, .2110,
-        .2802, .3589, .4443, 0.5368, 0.6255, 0.00, 0.00],  # CLI = 0.3
+        .2802, .3589, .4443, 0.5368, 0.6255, 0.00, 0.00, 0.00],  # CLI = 0.3
     [0.016, 0.020, .0294, .0478, .0678, .0893, .1118,
-        .1702, .2335, .3018, .3775, .4610, .5505, .6331],  # CLI = 0.4
+        .1702, .2335, .3018, .3775, .4610, .5505, .6331, 0.00],  # CLI = 0.4
     [0.00, 0.0324, .0486, .0671, .0875, .1094, .1326,
-        .1935, .2576, .3259, .3990, .4805, .5664, .6438],  # CLI = 0.5
-    [0.029, 0.043, 0.048, 0.049, .0524, .0684, .0868,
+        .1935, .2576, .3259, .3990, .4805, .5664, .6438, 0.00],  # CLI = 0.5
+    [0.00, 0.029, 0.043, 0.048, 0.049, .0524, .0684, .0868,
         .1074, .1298, .1537, .2169, .3512, .5025, .6605],  # CLI = 0.6
-    [.0510, .0743, .0891, .1074, .1281, .1509, .1753,
+    [0.00, .0510, .0743, .0891, .1074, .1281, .1509, .1753,
         .2407, .3083, .3775, .4496, .5265, .6065, .6826],  # CLI = 0.7
-    [.0670, .0973, .1114, .1290, .1494, .1723, .1972,
+    [0.00, .0670, .0973, .1114, .1290, .1494, .1723, .1972,
         .2646, .3345, .4047, .4772, .5532, .6307, .7092],  # CLI = 0.8
 ])
 CPEC = np.array(
@@ -394,24 +395,24 @@ def _biquad(T, i, xi, yi):
 # array length for CP_Angle_table and CT_Angle_table
 ang_arr_len = np.array([10, 6, 8, 8, 7, 10, 6])
 # array length for CP_CLi_table and CT_CLi_table
-cli_arr_len = np.array([12, 14, 14, 14, 14, 14])
+cli_arr_len = np.array([12, 14, 14, 15, 15, 15])
 # integrated design lift coefficient adjustment factor to power coefficient
 PF_CLI_arr = np.array([1.68, 1.405, 1.0, .655, .442, .255, .102])
 # integrated design lift coefficient adjustment factor to thrust coefficient
 TF_CLI_arr = np.array([1.22, 1.105, 1.0, .882, .792, .665, .540])
 num_blades_arr = np.array([2.0, 4.0, 6.0, 8.0])
 XPCLI = np.array([
     [4.26, 2.285, 1.780, 1.568, 1.452, 1.300, 1.220,
-        1.160, 1.110, 1.085, 1.054, 1.048, 0.000, 0.000],  # CL = 0.3
+        1.160, 1.110, 1.085, 1.054, 1.048, 0.000, 0.000, 0.0],  # CL = 0.3
     [2.0, 1.88, 1.652, 1.408, 1.292, 1.228, 1.188,
-        1.132, 1.105, 1.08, 1.058, 1.042, 1.029, 1.0220],  # CL = 0.4
+        1.132, 1.105, 1.08, 1.058, 1.042, 1.029, 1.0220, 0.0],  # CL = 0.4
     [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
-        1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000],  # CL = 0.5
-    [0.0, .40, .52, .551, .619, .712, .775,
+        1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 0.0],  # CL = 0.5
+    [0.0, 0.065, .40, .52, .551, .619, .712, .775,
         0.815, 0.845, 0.865, 0.891, 0.928, 0.958, 0.975],  # CL = 0.6
-    [0.00, .436, .545, .625, .682, .726, .755,
+    [0.00, 0.250, .436, .545, .625, .682, .726, .755,
         0.804, 0.835, 0.864, 0.889, 0.914, 0.935, 0.944],  # CL = 0.7
-    [0.00, .333, .436, .520, .585, .635, .670,
+    [0.00, 0.110, .333, .436, .520, .585, .635, .670,
         0.730, 0.770, 0.807, 0.835, 0.871, 0.897, 0.909],  # CL = 0.8
 ])
 XTCLI = np.array([
@@ -421,11 +422,11 @@ def _biquad(T, i, xi, yi):
         1.110, 1.089, 1.071, 1.060, 1.054, 1.051, 1.048],  # CL = 0.4
     [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
         1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000],  # CL = 0.5
-    [.000, .399, .694, .787, .831, .860, .881,
+    [.295, .399, .694, .787, .831, .860, .881,
         0.908, 0.926, 0.940, 0.945, 0.951, 0.958, 0.958],  # CL = 0.6
-    [.000, .251, .539, .654, .719, .760, .788,
+    [.166, .251, .539, .654, .719, .760, .788,
         0.831, 0.865, 0.885, 0.900, 0.910, 0.916, 0.916],  # CL = 0.7
-    [0.0, .1852, .442, .565, .635, .681, .716,
+    [0.042, .1852, .442, .565, .635, .681, .716,
         0.769, 0.809, 0.838, 0.855, 0.874, 0.881, 0.881],  # CL = 0.8
 ])
 adv_ratio_array2 = ([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])
@@ -550,7 +551,7 @@ class HamiltonStandard(om.ExplicitComponent):
     This is Hamilton Standard component rewritten from Fortran code. 
     The original documentation is available at 
     https://ntrs.nasa.gov/api/citations/19720010354/downloads/19720010354.pdf
-    It computes the thrust coefficient of a propeller.
+    It computes the thrust coefficient of a propeller blade.
     """
 
     def initialize(self):
@@ -577,6 +578,7 @@ def setup(self):
 
         self.add_output('thrust_coefficient', val=np.zeros(nn), units='unitless')
         self.add_output('ang_blade', val=np.zeros(nn), units='deg')
+        # Tip Compressibility loss factor
         self.add_output('comp_tip_loss_factor', val=np.zeros(nn), units='unitless')
 
         self.declare_partials('*', '*', method='fd', form='forward')
@@ -593,8 +595,8 @@ def compute(self, inputs, outputs):
             BLLL = np.zeros(7)
             PXCLI = np.zeros(7)
             XFFT = np.zeros(6)
-            CTG = np.zeros(6)
-            CTG1 = np.zeros(6)
+            CTG = np.zeros(11)
+            CTG1 = np.zeros(11)
             TXCLI = np.zeros(6)
             CTTT = np.zeros(4)
             XXXFT = np.zeros(4)
@@ -679,36 +681,39 @@ def compute(self, inputs, outputs):
                     CPE1 = CP_Eff*PBL*PF_CLI_arr[kdx]
                     CL_tab_idx = CL_tab_idx_begin
                     for kl in range(CL_tab_idx_begin, CL_tab_idx_end+1):
+                        CPE1X = CPE1
+                        if (CPE1 < CP_CLi_table[CL_tab_idx][0]):
+                            CPE1X = CP_CLi_table[CL_tab_idx][0]
                         cli_len = cli_arr_len[CL_tab_idx]
                         PXCLI[kl], run_flag = _unint(
-                            CP_CLi_table[CL_tab_idx][:cli_len], XPCLI[CL_tab_idx], CPE1)
+                            CP_CLi_table[CL_tab_idx][:cli_len], XPCLI[CL_tab_idx], CPE1X)
                         if (run_flag == 1):
                             ichck = ichck + 1
                         if (self.options[Settings.VERBOSITY] != Verbosity.DEBUG):
                             if (ichck <= 1):
                                 if (run_flag == 1):
                                     warnings.warn(
                                         f"Mach,VTMACH,J,power_coefficient,CP_Eff =: {inputs[Dynamic.Mission.MACH][i_node]},{inputs['tip_mach'][i_node]},{inputs['adv_ratio'][i_node]},{power_coefficient},{CP_Eff}")
-                                if (kl == 4 and CPE1 < 0.049):
+                                if (kl == 4 and CPE1 < 0.010):
                                     print(
                                         f"Extrapolated data is being used for CLI=.6--CPE1,PXCLI,L= , {CPE1},{PXCLI[kl]},{idx_blade}   Suggest inputting CLI=.5")
-                                if (kl == 5 and CPE1 < 0.0705):
+                                if (kl == 5 and CPE1 < 0.010):
                                     print(
                                         f"Extrapolated data is being used for CLI=.7--CPE1,PXCLI,L= , {CPE1},{PXCLI[kl]},{idx_blade}   Suggest inputting CLI=.5")
-                                if (kl == 6 and CPE1 < 0.0915):
+                                if (kl == 6 and CPE1 < 0.010):
                                     print(
                                         f"Extrapolated data is being used for CLI=.8--CPE1,PXCLI,L= , {CPE1},{PXCLI[kl]},{idx_blade}   Suggest inputting CLI=.5")
                         else:
                             if (run_flag == 1):
                                 warnings.warn(
                                     f"Mach,VTMACH,J,power_coefficient,CP_Eff =: {inputs[Dynamic.Mission.MACH][i_node]},{inputs['tip_mach'][i_node]},{inputs['adv_ratio'][i_node]},{power_coefficient},{CP_Eff}")
-                            if (kl == 4 and CPE1 < 0.049):
+                            if (kl == 4 and CPE1 < 0.010):
                                 print(
                                     f"Extrapolated data is being used for CLI=.6--CPE1,PXCLI,L= , {CPE1},{PXCLI[kl]},{idx_blade}   Suggest inputting CLI=.5")
-                            if (kl == 5 and CPE1 < 0.0705):
+                            if (kl == 5 and CPE1 < 0.010):
                                 print(
                                     f"Extrapolated data is being used for CLI=.7--CPE1,PXCLI,L= , {CPE1},{PXCLI[kl]},{idx_blade}   Suggest inputting CLI=.5")
-                            if (kl == 6 and CPE1 < 0.0915):
+                            if (kl == 6 and CPE1 < 0.010):
                                 print(
                                     f"Extrapolated data is being used for CLI=.8--CPE1,PXCLI,L= , {CPE1},{PXCLI[kl]},{idx_blade}   Suggest inputting CLI=.5")
                         NERPT = 1
@@ -746,19 +751,23 @@ def compute(self, inputs, outputs):
                 CTG[1] = .200
                 TFCLII, run_flag = _unint(
                     adv_ratio_array, TF_CLI_arr, inputs['adv_ratio'][i_node])
+                NCTG = 10
                 ifnd1 = 0
                 ifnd2 = 0
-                for il in range(5):
+                for il in range(NCTG):
                     ct = CTG[il]
                     CT_Eff = CTG[il]*AFCTE
                     TBL, run_flag = _unint(CTEC, BL_T_corr_table[idx_blade], CT_Eff)
                     # TBL = number of blades correction for thrust_coefficient
                     CTE1 = CT_Eff*TBL*TFCLII
                     CL_tab_idx = CL_tab_idx_begin
                     for kl in range(CL_tab_idx_begin, CL_tab_idx_end+1):
+                        CTE1X = CTE1
+                        if (CTE1 < CT_CLi_table[CL_tab_idx][0]):
+                            CTE1X = CT_CLi_table[CL_tab_idx][0]
                         cli_len = cli_arr_len[CL_tab_idx]
                         TXCLI[kl], run_flag = _unint(
-                            CT_CLi_table[CL_tab_idx][:cli_len], XTCLI[CL_tab_idx][:cli_len], CTE1)
+                            CT_CLi_table[CL_tab_idx][:cli_len], XTCLI[CL_tab_idx][:cli_len], CTE1X)
                         NERPT = 5
                         if (run_flag == 1):
                             # off lower bound only.

aviary/subsystems/propulsion/prop_performance.py

@@ -25,7 +25,7 @@ def setup(self):
         self.add_subsystem(
             name='sqa_comp',
             subsys=om.ExecComp(
-                'sqa = minimum(DiamNac**2/DiamProp**2, 0.32)',
+                'sqa = minimum(DiamNac**2/DiamProp**2, 0.50)',
                 DiamNac={'units': 'ft'},
                 DiamProp={'units': 'ft'},
                 sqa={'units': 'unitless'},
@@ -78,7 +78,7 @@ def setup(self):
         self.blockage_factor_interp.add_input(
             "sqa_array",
             0.0,
-            training_data=[0.00, 0.04, 0.08, 0.12, 0.16, 0.20, 0.24, 0.28, 0.32],
+            training_data=[0.00, 0.04, 0.08, 0.12, 0.16, 0.20, 0.24, 0.28, 0.32, 0.50],
             units="unitless",
             desc="square of DiamNac/DiamProp",
         )
@@ -105,7 +105,8 @@ def setup(self):
                  [0.964, 0.954, 0.943, 0.912, 0.876, 0.834, 0.786],
                  [0.955, 0.943, 0.928, 0.892, 0.848, 0.801, 0.751],
                  [0.948, 0.935, 0.917, 0.872, 0.820, 0.763, 0.706],
-                 [0.940, 0.924, 0.902, 0.848, 0.790, 0.726, 0.662]]
+                 [0.940, 0.924, 0.902, 0.848, 0.790, 0.726, 0.662],
+                 [0.904, 0.875, 0.835, 0.740, 0.655, 0.560, 0.464]]
             ),
         )
 

aviary/subsystems/propulsion/test/test_prop_performance.py

@@ -11,19 +11,21 @@
 from aviary.variable_info.options import get_option_defaults
 
 # Setting up truth values from GASP
-CT = np.array([0.27651, 0.20518, 0.13093, 0.10236,
-              0.10236, 0.19331, 0.10189, 0.10189, 0.18123])
-XFT = np.array([1.0, 1.0, 0.9976, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
-CTX = np.array([0.27651, 0.20518, 0.13062, 0.10236,
-               0.10236, 0.19331, 0.10189, 0.10189, 0.18123])
+CT = np.array([0.27651, 0.20518, 0.13093, 0.10236, 0.10236, 0.19331,
+               0.10189, 0.10189, 0.18123, 0.08523, 0.06463, 0.02800])
+XFT = np.array([1.0, 1.0, 0.9976, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
+CTX = np.array([0.27651, 0.20518, 0.13062, 0.10236, 0.10236, 0.19331,
+               0.10189, 0.10189, 0.18123, 0.08523, 0.06463, 0.02800])
 three_quart_blade_angle = np.array(
-    [25.17, 29.67, 44.23, 31.94, 31.94, 17.44, 33.43, 33.43, 20.08])
-thrust = np.array([4634.8, 3415.9, 841.5, 1474.3, 1400.6,
-                  3923.5, 1467.6, 1394.2, 3678.3])
-prop_eff = np.array([0.0, 0.7235, 0.892, 0.9059, 0.9059, 0.5075, 0.9017, 0.9017, 0.4758])
-install_loss = np.array([0.0133, 0.02, 0.034, 0.0, 0.05, 0.05, 0.0, 0.05, 0.05])
-install_eff = np.array([0.0, 0.709, 0.8617, 0.9059, 0.8606,
-                       0.4821, 0.9017, 0.8566, 0.452])
+    [25.17, 29.67, 44.23, 31.94, 31.94, 17.44, 33.43, 33.43, 20.08, 30.28, 29.50, 28.10])
+thrust = np.array([4634.8, 3415.9, 841.5, 1474.3, 1400.6, 3923.5,
+                   1467.6, 1394.2, 3678.3, 1210.4, 917.8, 397.7])
+prop_eff = np.array([0.00078, 0.72352, 0.89202, 0.90586, 0.90586, 0.50750,
+                     0.90172, 0.90172, 0.47579, 0.83809, 0.76259, 0.49565])
+install_loss = np.array([0.0133, 0.02, 0.034, 0.0, 0.05, 0.05,
+                         0.0, 0.05, 0.05, 0.0140, 0.0140, 0.0140])
+install_eff = np.array([0.00077, 0.70904, 0.86171, 0.90586, 0.86056, 0.48213,
+                        0.90172, 0.85664, 0.45200, 0.82635, 0.75190, 0.48871])
 
 
 class PropPerformanceTest(unittest.TestCase):
@@ -169,6 +171,34 @@ def test_case_6_7_8(self):
             minimum_step=1e-12, abs_err_tol=5.0E-4, rel_err_tol=5.0E-5, excludes=["*atmosphere*"])
         assert_check_partials(partial_data, atol=1e-4, rtol=1e-4)
 
+    def test_case_9_10_11(self):
+        # Case 9, 10, 11, to test CLI > 0.5
+        prob = self.prob
+        prob.set_val(Aircraft.Engine.PROPELLER_DIAMETER, 12.0, units="ft")
+        prob.set_val(Aircraft.Nacelle.AVG_DIAMETER, 2.4, units='ft')
+        prob.set_val(Aircraft.Engine.PROPELLER_ACTIVITY_FACTOR, 150.0, units="unitless")
+        prob.set_val(Aircraft.Engine.PROPELLER_INTEGRATED_LIFT_COEFFICENT,
+                     0.65, units="unitless")
+        prob.set_val(Dynamic.Mission.ALTITUDE, [10000.0, 10000.0, 10000.0], units="ft")
+        prob.set_val(Dynamic.Mission.VELOCITY, [200.0, 200.0, 200.0], units="knot")
+        prob.set_val(Dynamic.Mission.PROPELLER_TIP_SPEED,
+                     [750.0, 750.0, 750.0], units="ft/s")
+        prob.set_val(Dynamic.Mission.SHAFT_POWER, [900.0, 750.0, 500.0], units="hp")
+
+        prob.run_model()
+        self.compare_results(case_idx_begin=9, case_idx_end=11)
+
+        partial_data = prob.check_partials(
+            out_stream=None, compact_print=True, show_only_incorrect=True, form='central', method="fd",
+            minimum_step=1e-12, abs_err_tol=5.0E-4, rel_err_tol=5.0E-5, excludes=["*atmosphere*"])
+        # remove partial derivative of 'comp_tip_loss_factor' with respect to
+        # 'aircraft:engine:propeller_integrated_lift_coefficient' from assert_check_partials
+        partial_data_hs = partial_data['pp.hamilton_standard']
+        key_pair = ('comp_tip_loss_factor',
+                    'aircraft:engine:propeller_integrated_lift_coefficient')
+        del partial_data_hs[key_pair]
+        assert_check_partials(partial_data, atol=1.5e-3, rtol=1e-4)
+
 
 if __name__ == "__main__":
     unittest.main()

aviary/variable_info/variable_meta_data.py

@@ -6173,12 +6173,12 @@
 add_meta_data(
     Dynamic.Mission.PROPELLER_TIP_SPEED,
     meta_data=_MetaData,
-    historical_name={"GASP": 'INGASP.TSPDMX',
+    historical_name={"GASP": None,
                      "FLOPS": None,
                      "LEAPS1": None
                      },
     units='ft/s',
-    desc='maximum allowable propeller tip speed',
+    desc='propeller tip speed',
     default_value=500.0,
 )