Remove vectorize and expand natural convection range

docs/refs.bib

@@ -22,6 +22,15 @@ @techreport{cigre207
   year        = {2002}
 }
 
+@book{incropera2007,
+  address = {New York City, New York},
+  author = {Incropera, Frank P. and DeWitt, David P. and Bergman, Theodore L. and Lavine, Adrienne S.},
+  edition = {6th Edition},
+  publisher = {John Wiley & Sons, Inc.},
+  title = {Fundamentals of Heat and Mass Transfer},
+  year = 2007
+}
+
 @ARTICLE{ieee.acsr.taskforce,
   author={Rathbun, L.S. and Douglass, D.A. and Kirkpatrick, L.A.},
   journal={IEEE Transactions on Power Apparatus and Systems},

linerate/equations/cigre207/convective_cooling.py

@@ -1,7 +1,6 @@
 import warnings
 
 import numpy as np
-from numba import vectorize
 
 from ...units import (
     Celsius,
@@ -14,6 +13,7 @@
     WattPerMeterPerKelvin,
 )
 
+
 # Physical quantities
 #####################
 
@@ -82,7 +82,7 @@ def compute_kinematic_viscosity_of_air(film_temperature: Celsius) -> KilogramPer
 
 
 def compute_prandtl_number(
-    film_temperature: Celsius,
+        film_temperature: Celsius,
 ) -> Unitless:
     r"""Compute the Prandtl number.
 
@@ -106,10 +106,10 @@ def compute_prandtl_number(
 
 
 def compute_reynolds_number(
-    wind_speed: MeterPerSecond,
-    conductor_diameter: Meter,
-    kinematic_viscosity_of_air: SquareMeterPerSecond,
-    relative_air_density: Unitless,
+        wind_speed: MeterPerSecond,
+        conductor_diameter: Meter,
+        kinematic_viscosity_of_air: SquareMeterPerSecond,
+        relative_air_density: Unitless,
 ) -> Unitless:
     r"""Compute the Reynolds number using the conductor diameter as characteristic length scale.
 
@@ -153,30 +153,24 @@ def _check_perpendicular_flow_nusseltnumber_out_of_bounds(reynolds_number):
         warnings.warn("Reynolds number is out of bounds", stacklevel=5)
 
 
-@vectorize(nopython=True)
 def _compute_perpendicular_flow_nusseltnumber(
-    reynolds_number: Unitless,
-    conductor_roughness: Meter,
+        reynolds_number: Unitless,
+        conductor_roughness: Meter,
 ) -> Unitless:
     # From table on page 6 in Cigre207
     Re = reynolds_number
     Rs = conductor_roughness
-
-    if Re < 100:
-        B, n = 0, 0
-    elif Re < 2.65e3:
-        B, n = 0.641, 0.471
-    elif Rs <= 0.05:
-        B, n = 0.178, 0.633
-    else:
-        B, n = 0.048, 0.800
-
-    return B * Re**n  # type: ignore
+    conditions = [Re < 100, Re < 2.65e3, Rs <= 0.05]
+    B_choices = [0, 0.641, 0.178]
+    n_choices = [0, 0.471, 0.633]
+    B = np.select(conditions, B_choices, default=0.048)
+    n = np.select(conditions, n_choices, default=0.800)
+    return B * Re ** n  # type: ignore
 
 
 def compute_perpendicular_flow_nusseltnumber(
-    reynolds_number: Unitless,
-    conductor_roughness: Meter,
+        reynolds_number: Unitless,
+        conductor_roughness: Meter,
 ) -> Unitless:
     r"""Compute the Nusselt number for perpendicular flow.
 
@@ -202,7 +196,7 @@ def compute_perpendicular_flow_nusseltnumber(
 
 
 def compute_low_wind_speed_nusseltnumber(
-    perpendicular_flow_nusselt_number: Unitless,
+        perpendicular_flow_nusselt_number: Unitless,
 ) -> Unitless:
     r"""Compute the corrected Nusselt number for low wind speed.
 
@@ -219,27 +213,9 @@ def compute_low_wind_speed_nusseltnumber(
     return 0.55 * perpendicular_flow_nusselt_number
 
 
-@vectorize(nopython=True)
-def _correct_wind_direction_effect_on_nusselt_number(
-    perpendicular_flow_nusselt_number: Unitless,
-    angle_of_attack: Radian,
-) -> Unitless:
-    delta = angle_of_attack
-    Nu_90 = perpendicular_flow_nusselt_number
-
-    sin_delta = np.sin(delta)
-    # Equation (14) on page 7 of Cigre207
-    if delta <= np.radians(24):
-        correction_factor = 0.42 + 0.68 * (sin_delta**1.08)
-    else:
-        correction_factor = 0.42 + 0.58 * (sin_delta**0.90)
-
-    return correction_factor * Nu_90
-
-
 def correct_wind_direction_effect_on_nusselt_number(
-    perpendicular_flow_nusselt_number: Unitless,
-    angle_of_attack: Radian,
+        perpendicular_flow_nusselt_number: Unitless,
+        angle_of_attack: Radian,
 ) -> Unitless:
     r"""Correct the Nusselt number for the wind's angle-of-attack.
 
@@ -260,54 +236,54 @@ def correct_wind_direction_effect_on_nusselt_number(
     Union[float, float64, ndarray[Any, dtype[float64]]]
         :math:`\text{Nu}_\delta`. The Nusselt number for the given wind angle-of-attack.
     """
-    return _correct_wind_direction_effect_on_nusselt_number(
-        perpendicular_flow_nusselt_number, angle_of_attack
-    )
+    delta = angle_of_attack
+    Nu_90 = perpendicular_flow_nusselt_number
+    sin_delta = np.sin(delta)
+    # Equation (14) on page 7 of Cigre207
+    return (np.where(delta <= np.radians(24),
+                     0.42 + 0.68 * (sin_delta ** 1.08),
+                     0.42 + 0.58 * (sin_delta ** 0.90))
+            * Nu_90)
 
 
 ## Natural convection computations (no wind):
 #############################################
 
 
 def _check_horizontal_natural_nusselt_number(
-    grashof_number: Unitless, prandtl_number: Unitless
+        grashof_number: Unitless, prandtl_number: Unitless
 ) -> None:
     GrPr = grashof_number * prandtl_number
     if np.any(GrPr < 0):
         raise ValueError("GrPr cannot be negative.")
-    elif np.any(GrPr > 1e6):
-        raise ValueError("GrPr out of bounds: Must be < 10^6.")
+    elif np.any(GrPr >= 1e12):
+        raise ValueError("GrPr out of bounds: Must be < 10^12.")
 
 
-@vectorize(nopython=True)
 def _compute_horizontal_natural_nusselt_number(
-    grashof_number: Unitless,
-    prandtl_number: Unitless,
+        grashof_number: Unitless,
+        prandtl_number: Unitless,
 ) -> Unitless:
     GrPr = grashof_number * prandtl_number
-
-    if GrPr < 1e2:
-        # Outside table range, should we use 0??
-        return 0
-    elif GrPr < 1e4:
-        return 0.850 * GrPr**0.188
-    elif GrPr < 1e6:
-        return 0.480 * GrPr**0.250
-    else:
-        # Outside table range, what should we do here?
-        return 0.125 * GrPr**0.333
+    # CIGRE 207 only allows GrPr in the range 1e2-1e6.
+    # In Incropera (2006), Table 9.1, the same values appear, but the table covers a wider range
+    # from 1e-10 to 1e12, which we use here
+    conditions = [GrPr < 1e-2, GrPr < 1e2, GrPr < 1e4, GrPr < 1e7, GrPr < 1e12]
+    n_choices = [0.058, 0.148, 0.188, 0.250, 0.333]
+    C_choices = [0.0675, 1.02, 0.850, 0.480, 0.125]
+    C = np.select(conditions, C_choices, default=np.nan)
+    n = np.select(conditions, n_choices, default=np.nan)
+    return C * GrPr ** n  # type: ignore
 
 
 def compute_horizontal_natural_nusselt_number(
-    grashof_number: Unitless,
-    prandtl_number: Unitless,
+        grashof_number: Unitless,
+        prandtl_number: Unitless,
 ) -> Unitless:
     r"""The Nusselt number for natural (passive) convection on a horizontal conductor.
 
     Equation (16) and Table II on page 7 of :cite:p:`cigre207`.
-
-    The coefficient table is modified slightly so coefficients with
-    :math:`\text{Gr}\text{Pr} < 0.1` leads to :math:`\text{Nu} = 0`.
+    We expand the allowable range by using Table 9.1 of :cite:p:`incropera2007`.
 
     Parameters
     ----------
@@ -329,10 +305,10 @@ def compute_horizontal_natural_nusselt_number(
 
 
 def compute_nusselt_number(
-    forced_convection_nusselt_number: Unitless,
-    natural_nusselt_number: Unitless,
-    low_wind_nusselt_number: Unitless,
-    wind_speed: MeterPerSecond,
+        forced_convection_nusselt_number: Unitless,
+        natural_nusselt_number: Unitless,
+        low_wind_nusselt_number: Unitless,
+        wind_speed: MeterPerSecond,
 ) -> Unitless:
     r"""Compute the nusselt number.