Initial implementation of CIGRE 207

linerate/equations/cigre207/__init__.py

@@ -0,0 +1,5 @@
+"""
+This submodule contains implementations of equations listed in :cite:p:`cigre207`.
+"""
+
+from . import convective_cooling, solar_heating  # noqa

linerate/equations/cigre207/ac_resistance.py

@@ -0,0 +1,45 @@
+from typing import Union
+
+from linerate.units import OhmPerMeter, Ampere, SquareMeter, Unitless, SquareMeterPerAmpere
+
+
+def correct_resistance_acsr_magnetic_core_loss_simple(
+    ac_resistance: OhmPerMeter,
+    current: Ampere,
+    aluminium_cross_section_area: SquareMeter,
+    constant_magnetic_effect: Union[Unitless, None],
+    current_density_proportional_magnetic_effect: Union[SquareMeterPerAmpere, None],
+    max_relative_increase: Unitless,
+) -> OhmPerMeter:
+    r"""
+    Return resistance with constant correction for magnetic effects, using simple method from
+    Cigre 207, see section 2.1.1.
+
+    Parameters
+    ----------
+    ac_resistance:
+        :math:`R~\left[\Omega\right]`. The AC resistance of the conductor.
+    current:
+        :math:`I~\left[\text{A}\right]`. The current going through the conductor.
+    aluminium_cross_section_area:
+        :math:`A_{\text{Al}}~\left[\text{m}^2\right]`. The cross sectional area of the aluminium
+        strands in the conductor.
+    constant_magnetic_effect:
+        :math:`b`. The constant magnetic effect, most likely equal to 1. If ``None``, then no
+        correction is used (useful for non-ACSR cables).
+    current_density_proportional_magnetic_effect:
+        :math:`m`. The current density proportional magnetic effect. If ``None``, then it is
+        assumed equal to 0.
+    max_relative_increase:
+        :math:`c_\text{max}`. Saturation point of the relative increase in conductor resistance.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`R_\text{corrected}~\left[\Omega\right]`. The resistance of the conductor after
+        taking steel core magnetization effects into account.
+    """
+    return 1.0123 * ac_resistance
+
+
+# TODO: Implement section 2.1.2?

linerate/equations/cigre207/convective_cooling.py

@@ -0,0 +1,331 @@
+import warnings
+
+import numpy as np
+from numba import vectorize
+
+from ...units import (
+    Celsius,
+    KilogramPerCubeMeter,
+    Meter,
+    MeterPerSecond,
+    Radian,
+    Unitless,
+    WattPerMeterPerKelvin,
+)
+
+
+# Physical quantities
+#####################
+
+
+def compute_thermal_conductivity_of_air(film_temperature: Celsius) -> WattPerMeterPerKelvin:
+    r"""Approximation of the thermal conductivity of air.
+
+    On page 5 of :cite:p:`cigre207`.
+
+    Parameters
+    ----------
+    film_temperature:
+        :math:`T_f = 0.5 (T_s + T_a)~\left[^\circ\text{C}\right]`. The temperature of the
+        thin air-film surrounding the conductor. Equal to the average of the ambient air
+        temperature and the conductor sufrace temperature.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`\lambda_f~\left[\text{W}~\text{m}^{-1}~\text{K}^{-1}\right]`. The thermal
+        conductivity of air at the given temperature.
+    """
+    T_f = film_temperature
+    return 2.42e-2 + 7.2e-5 * T_f
+
+
+def compute_relative_air_density(
+    height_above_sea_level: Meter
+) -> Unitless:
+    r"""Approximation of the relative density of air at a given altitude, relative to density at sea level.
+
+    Equation on page 6 of :cite:p:`cigre207`.
+
+    Parameters
+    ----------
+    height_above_sea_level:
+        :math:`y~\left[\text{m}\right]`. The conductor's altitude.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`\rho_r`. The relative mass density of air.
+    """
+    y = height_above_sea_level
+    return np.exp(-1.16e-4*y)
+
+
+def compute_kinematic_viscosity_of_air(film_temperature: Celsius) -> KilogramPerCubeMeter:
+    r"""Approximation of the kinematic viscosity of air at a given temperature.
+
+    Equation on page 5 of :cite:p:`cigre207`.
+
+    Parameters
+    ----------
+    film_temperature:
+        :math:`T_f = 0.5 (T_s + T_a)~\left[^\circ\text{C}\right]`. The temperature of the
+        thin air-film surrounding the conductor. Equal to the average of the ambient air
+        temperature and the conductor sufrace temperature.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`\nu_f~\left[\text{m}^2~\text{s}^{-1}\right]`. The kinematic viscosity of air.
+    """
+    T_f = film_temperature
+    return 1.32e-5 + 9.5e-8*T_f
+
+
+
+def compute_prandtl_number(
+    film_temperature: Celsius,
+) -> Unitless:
+    r"""Compute the Prandtl number.
+
+    Defined on page 5 of :cite:p:`cigre207`.
+
+    The Prandtl number measures the ratio between viscosity and thermal diffusivity for a fluid.
+
+    Parameters
+    ----------
+    film_temperature:
+        :math:`T_f = 0.5 (T_s + T_a)~\left[^\circ\text{C}\right]`. The temperature of the
+        thin air-film surrounding the conductor. Equal to the average of the ambient air
+        temperature and the conductor sufrace temperature.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`\text{Pr}`. The Prandtl number.
+    """
+    return 0.715 - 2.5e-4*film_temperature
+
+
+## Nusselt number calculation
+#############################
+
+
+def _check_perpendicular_flow_nusseltnumber_out_of_bounds(reynolds_number):
+    Re = reynolds_number
+    if np.any(Re < 0):
+        raise ValueError("Reynolds number cannot be negative.")
+
+    if np.any(np.logical_or(Re < 100, Re > 5e4)):
+        warnings.warn("Reynolds number is out of bounds", stacklevel=5)
+
+
+@vectorize(nopython=True)
+def _compute_perpendicular_flow_nusseltnumber(
+    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
+
+
+def compute_perpendicular_flow_nusseltnumber(
+    reynolds_number: Unitless,
+    conductor_roughness: Meter,
+) -> Unitless:
+    r"""Compute the Nusselt number for perpendicular flow.
+
+    The Nusselt number is the ratio of conductive heat transfer to convective heat transfer.
+
+    Parameters
+    ----------
+    reynolds_number:
+        :math:`\text{Re}`. The Reynolds number.
+    conductor_roughness:
+        :math:`\text{Rs}`. The roughness number
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`\text{Nu}_{90}`. The perpendicular flow Nusselt number.
+    """
+    _check_perpendicular_flow_nusseltnumber_out_of_bounds(reynolds_number)
+    return _compute_perpendicular_flow_nusseltnumber(
+        reynolds_number,
+        conductor_roughness,
+    )
+
+
+def compute_low_wind_speed_nusseltnumber(
+    perpendicular_flow_nusselt_number: Unitless,
+) -> Unitless:
+    r"""Compute the corrected Nusselt number for low wind speed.
+
+    Parameters
+    ----------
+    perpendicular_flow_nusselt_number:
+        :math:`\text{Nu}_{90}`. The perpendicular flow Nusselt number.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`\text{Nu}_{90}`. The corrected Nusselt number for low wind speed.
+    """
+    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,
+) -> Unitless:
+    r"""Correct the Nusselt number for the wind's angle-of-attack.
+
+    Equation (14) on page 7 of :cite:p:`cigre207`.
+
+    The perpendicular flow nusselt number is denoted as :math:`\text{Nu}_\delta` in
+    :cite:p:`cigre207` since the wind's angle of attack is :math:`\delta`.
+
+    Parameters
+    ----------
+    perpendicular_flow_nusselt_number:
+        :math:`\text{Nu}_{90}`. The perpendicular flow Nusselt number.
+    angle_of_attack:
+        :math:`\delta~\left[\text{radian}\right]`. The wind angle-of-attack.
+
+    Returns
+    -------
+    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
+    )
+
+
+## Natural convection computations (no wind):
+#############################################
+
+
+def _check_horizontal_natural_nusselt_number(
+    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.")
+
+
+@vectorize(nopython=True)
+def _compute_horizontal_natural_nusselt_number(
+    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:
+        warnings.warn("GrPr out of bounds: Must be < 10^6.")
+        # Outside table range, what should we do here?
+        return 0.125 * GrPr**0.333
+
+
+def compute_horizontal_natural_nusselt_number(
+    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`.
+
+    Parameters
+    ----------
+    grashof_number:
+        :math:`\text{Gr}`. The Grashof number.
+    prandtl_number:
+        :math:`\text{Pr}`. The Prandtl number.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`\text{Nu}_0`. The natural convection nusselt number assuming horizontal conductor.
+    """
+    _check_horizontal_natural_nusselt_number(grashof_number, prandtl_number)
+    return _compute_horizontal_natural_nusselt_number(
+        grashof_number,
+        prandtl_number,
+    )
+
+
+def compute_nusselt_number(
+    forced_convection_nusselt_number: Unitless,
+    natural_nusselt_number: Unitless,
+    low_wind_nusselt_number: Unitless,
+    wind_speed: MeterPerSecond
+) -> Unitless:
+    r"""Compute the nusselt number.
+
+    Described on page 7 of :cite:p:`cigre207`.
+
+    Parameters
+    ----------
+    forced_convection_nusselt_number:
+        :math:`\text{Nu}_\delta`. The Nusselt number for the given wind angle-of-attack.
+    natural_nusselt_number:
+        :math:`\text{Nu}_\delta`. The natural convection nusselt number for horizontal conductor.
+    low_wind_nusselt_number:
+        :math:`\text{Nu}_cor`. Corrected Nusselt number for low wind.
+    wind_speed:
+        :math:`v~\left[\text{m}~\text{s}^{-1}\right]`. The wind speed.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`Nu`. The nusselt number.
+    """
+    max_nusselt = np.maximum(forced_convection_nusselt_number, natural_nusselt_number)
+    if wind_speed < 0.5:
+        return np.maximum(max_nusselt, low_wind_nusselt_number)
+    else:
+        return max_nusselt