From 78645edaeb89af6382593f05a863ae8e383fe211 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Tue, 28 May 2024 14:54:52 +0200 Subject: [PATCH 01/18] Initial implementation of CIGRE 207 --- linerate/equations/cigre207/__init__.py | 5 + linerate/equations/cigre207/ac_resistance.py | 45 +++ .../equations/cigre207/convective_cooling.py | 331 ++++++++++++++++++ linerate/equations/cigre207/solar_heating.py | 216 ++++++++++++ .../equations/cigre601/convective_cooling.py | 301 +++------------- linerate/equations/convective_cooling.py | 40 +++ linerate/equations/dimensionless.py | 166 +++++++++ .../equations/ieee738/convective_cooling.py | 45 --- linerate/model.py | 14 +- linerate/models/Cigre207.md | 30 ++ linerate/models/cigre207.py | 117 +++++++ 11 files changed, 1004 insertions(+), 306 deletions(-) create mode 100644 linerate/equations/cigre207/__init__.py create mode 100644 linerate/equations/cigre207/ac_resistance.py create mode 100644 linerate/equations/cigre207/convective_cooling.py create mode 100644 linerate/equations/cigre207/solar_heating.py create mode 100644 linerate/equations/convective_cooling.py create mode 100644 linerate/equations/dimensionless.py create mode 100644 linerate/models/Cigre207.md create mode 100644 linerate/models/cigre207.py diff --git a/linerate/equations/cigre207/__init__.py b/linerate/equations/cigre207/__init__.py new file mode 100644 index 0000000..af02088 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/linerate/equations/cigre207/ac_resistance.py b/linerate/equations/cigre207/ac_resistance.py new file mode 100644 index 0000000..86fa2ff --- /dev/null +++ b/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? \ No newline at end of file diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py new file mode 100644 index 0000000..e3b12f7 --- /dev/null +++ b/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 diff --git a/linerate/equations/cigre207/solar_heating.py b/linerate/equations/cigre207/solar_heating.py new file mode 100644 index 0000000..5f0cf88 --- /dev/null +++ b/linerate/equations/cigre207/solar_heating.py @@ -0,0 +1,216 @@ +import numpy as np +from numpy import pi + +from ...units import Meter, Unitless, WattPerMeter, WattPerSquareMeter + + +def compute_direct_solar_radiation( + sin_solar_altitude: Unitless, + clearness_ratio: Unitless, + height_above_sea_level: Meter, +) -> WattPerSquareMeter: + r"""Compute the direct solar radiation. + + Equation (10-11) on page 19 of :cite:p:`cigre601`. Equation (10) states that the direct solar + radiation on a surface normal to the solar beam at sea level, :math:`I_{B(0)}`, is given by + + .. math:: + + N_s \frac{1280 \sin(H_s)}{\sin(H_s) + 0.314}, + + where :math:`N_s` is the clearness ratio which is used to adjust the amount of radiation + compared to what goes through a standard Indian atmosphere, and :math:`H_s` is the solar + altitude. + + While the solar radiation model is based on :cite:p:`sharma1965interrelationships` and + therefore have parameters estimated for an Indian atmosphere, it gives comparable results to + the solar radiation model in the IEEE standard :cite:p:`ieee738`. It is therefore reasonable to + assume that the parameters work in other climates as well. + + Parameters + ---------- + sin_solar_altitude: + :math:`\sin\left(H_s\right)`. The sine of the solar altitude. + clearness_ratio: + :math:`N_s`. The clearness ratio (or clearness number in + :cite:p:`sharma1965interrelationships,cigre207`). + height_above_sea_level: + :math:`y~\left[\text{m}\right]`. The conductor's altitude. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. + + Note + ---- + The 1280 originates and 0.314 in the above equation originates from + :cite:p:`sharma1965interrelationships`, which is cited in :cite:p:`morgan1982thermal` (which is + listed as the reference in :cite:p:`cigre601`). In :cite:p:`sharma1965interrelationships` the + empirical relationship + + .. math:: + + I_{B(0)} = \frac{1.842 \sin(H_s)}{\sin(H_s) + 0.3135}~\text{Ly}~\text{min}^{-1} + + is introduced, and by converting from Langley per minute to :math:`\text{W}~\text{m}^{-2}`, we + obtain + + .. math:: + + I_{B(0)} = N_s \frac{1284.488 \sin(H_s)}{\sin(H_s) + 0.3135}~\text{W}~\text{m}^{-2}, + + which is equal to the equation we use (with three significant digits). + """ + sin_H_s = sin_solar_altitude + N_s = clearness_ratio + y = height_above_sea_level + + I_B_0 = N_s * 1280 * sin_H_s / (sin_H_s + 0.314) + # Equation 19 says that + # I_B = I_B_0 * (1 + 1.4e-4 * y * (1367/I_B_0 - 1)) + # However, if I_B_0 = 0, this will divide by 0. To return NaN-values if and only + # if the input is NaN, we therefore reformulate it + scaled_y = 1.4e-4 * y + I_B = I_B_0 * (1 - scaled_y) + 1367 * scaled_y + + return np.where(sin_H_s >= 0, I_B, 0 * I_B) + + +def compute_diffuse_sky_radiation( + direct_solar_radiation: WattPerSquareMeter, + sin_solar_altitude: Unitless, +) -> WattPerSquareMeter: + r"""Compute the diffuse radiation (light scattered in the atmosphere). + + Equation (13) on page 20 of :cite:p:`cigre601`. + + This equation differ from :cite:p:`cigre207`, however the difference is small, and the + diffuse radiation is a small contributor to the overall solar radiation, so the total + discrepancy between the models is small. + + Parameters + ---------- + direct_solar_radiation: + :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. + sin_solar_altitude: + :math:`\sin\left(H_s\right)`. The sine of the solar altitude. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I_d~\left[\text{W}~\text{m}^{-2}\right]`.The diffuse solar radiation. + """ + sin_H_s = sin_solar_altitude + I_B = direct_solar_radiation + return np.maximum(0, (430.5 - 0.3288 * I_B)) * np.maximum(0, sin_H_s) + + +def compute_global_radiation_intensity( + direct_solar_radiation: WattPerSquareMeter, + diffuse_sky_radiation: WattPerSquareMeter, + albedo: Unitless, + sin_angle_of_sun_on_line: Unitless, + sin_solar_altitude: Unitless, +) -> WattPerSquareMeter: + r"""Compute the global radiation intensity experienced by the conductor. + + Equation (9) on page 18 of :cite:p:`cigre601` state that the global radiation intensity, + :math:`I_T`, is given by + + .. math:: + + I_T = + I_B \left(\sin(\eta) + 0.5 F \pi \sin(H_s)\right) + + I_d \left(1 + 0.5 F \pi\right), + + where :math:`\eta` is the incidence angle of the sun on the line, :math:`H_s` is the solar + altitude and :math:`F` is the ground albedo (amount of radiation diffusely reflected from the + ground). The factor :math:`0.5 \pi` is due the assumption that the ground reflects light + diffusely and uniformly in all directions, so the reflected energy is always directed normally + to the line. In CIGRE207, it is also assumed that the diffuse radiation is uniformly directed, + which leads to :math:`I_d (0.5 \pi + 0.5 F \pi)` instead of :math:`I_d (1 + 0.5 F \pi)` + + Parameters + ---------- + direct_solar_radiation: + :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. + diffuse_sky_radiation: + :math:`I_d~\left[\text{W}~\text{m}^{-2}\right]`.The diffuse solar radiation. + albedo: + :math:`F`. The ground albedo. + sin_angle_of_sun_on_line: + :math:`\sin\left(\eta\right)`. The sine of the angle of the sun on the line. + sin_solar_altitude: + :math:`\sin\left(H_s\right)`. The sine of the solar altitude. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I_T~\left[\text{W}~\text{m}^{-2}\right]`. The global radiation intensity. + + Note + ---- + The following values are given for the albedo in :cite:p:`cigre601`: + + .. list-table:: + :widths: 50 50 + :header-rows: 1 + + * - Ground + - Albedo + * - Water (:math:`H_s > 30^\circ`) + - 0.05 + * - Forest + - 0.1 + * - Urban areas + - 0.15 + * - Soil, grass and crops + - 0.2 + * - Sand + - 0.3 + * - Ice + - 0.4-0.6 + * - Snow + - 0.6-0.8 + """ + I_B = direct_solar_radiation + I_d = diffuse_sky_radiation + F = albedo + sin_H_s = sin_solar_altitude + sin_eta = sin_angle_of_sun_on_line + F_pi_half = 0.5 * pi * F + + return I_B * (sin_eta + F_pi_half * sin_H_s) + I_d * (1 + F_pi_half) # type: ignore + + +def compute_solar_heating( + absorptivity: Unitless, + global_radiation_intensity: WattPerSquareMeter, + conductor_diameter: Meter, +) -> WattPerMeter: + r"""Compute the solar heating experienced by the conductor. + + Equation (8) on page 18 of :cite:p:`cigre601`. + + Parameters + ---------- + absorptivity: + :math:`\alpha_s`. Material constant. According to :cite:p:`cigre601`, it starts at + approximately 0.2 for new cables and reaches a constant value of approximately 0.9 + after about one year. + global_radiation_intensity: + :math:`I_T~\left[\text{W}~\text{m}^{-2}\right]`.The global radiation intensity. + conductor_diameter: + :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_S~\left[\text{W}~\text{m}^{-1}\right]`. The solar heating of the conductor + """ + alpha_s = absorptivity + I_T = global_radiation_intensity + D = conductor_diameter + + return alpha_s * I_T * D diff --git a/linerate/equations/cigre601/convective_cooling.py b/linerate/equations/cigre601/convective_cooling.py index b4e70bf..c02ff29 100644 --- a/linerate/equations/cigre601/convective_cooling.py +++ b/linerate/equations/cigre601/convective_cooling.py @@ -1,17 +1,13 @@ import warnings -from textwrap import dedent import numpy as np from numba import vectorize from ...units import ( Celsius, - JoulePerKilogramPerKelvin, KilogramPerCubeMeter, KilogramPerMeterPerSecond, Meter, - MeterPerSecond, - MeterPerSquareSecond, Radian, SquareMeterPerSecond, Unitless, @@ -19,15 +15,16 @@ WattPerMeterPerKelvin, ) + # Physical quantities ##################### def compute_temperature_gradient( - total_heat_gain: WattPerMeter, - conductor_thermal_conductivity: WattPerMeterPerKelvin, - core_diameter: Meter, - conductor_diameter: Meter, + total_heat_gain: WattPerMeter, + conductor_thermal_conductivity: WattPerMeterPerKelvin, + core_diameter: Meter, + conductor_diameter: Meter, ) -> Celsius: r"""Compute the difference between the core and surface temperature. @@ -67,8 +64,8 @@ def compute_temperature_gradient( if D_1 == 0: # TODO: Maybe lower tolerance? return 0.5 * tmp else: - D_1_sq = D_1**2 - delta_D_sq = D**2 - D_1_sq + D_1_sq = D_1 ** 2 + delta_D_sq = D ** 2 - D_1_sq return tmp * (0.5 - (D_1_sq / delta_D_sq) * np.log(D / D_1)) @@ -93,11 +90,11 @@ def compute_thermal_conductivity_of_air(film_temperature: Celsius) -> WattPerMet conductivity of air at the given temperature. """ T_f = film_temperature - return 2.368e-2 + 7.23e-5 * T_f - 2.763e-8 * (T_f**2) + return 2.368e-2 + 7.23e-5 * T_f - 2.763e-8 * (T_f ** 2) def compute_air_density( - film_temperature: Celsius, height_above_sea_level: Meter + film_temperature: Celsius, height_above_sea_level: Meter ) -> KilogramPerCubeMeter: r"""Approximation of the density of air at a given temperature and altitude. @@ -119,7 +116,7 @@ def compute_air_density( """ T_f = film_temperature y = height_above_sea_level - return (1.293 - 1.525e-4 * y + 6.379e-9 * (y**2)) / (1 + 0.00367 * T_f) + return (1.293 - 1.525e-4 * y + 6.379e-9 * (y ** 2)) / (1 + 0.00367 * T_f) def compute_dynamic_viscosity_of_air(film_temperature: Celsius) -> KilogramPerMeterPerSecond: @@ -141,11 +138,11 @@ def compute_dynamic_viscosity_of_air(film_temperature: Celsius) -> KilogramPerMe of air. """ T_f = film_temperature - return 17.239e-6 + 4.635e-8 * T_f - 2.03e-11 * (T_f**2) + return 17.239e-6 + 4.635e-8 * T_f - 2.03e-11 * (T_f ** 2) def compute_kinematic_viscosity_of_air( - dynamic_viscosity_of_air: KilogramPerMeterPerSecond, air_density: KilogramPerCubeMeter + dynamic_viscosity_of_air: KilogramPerMeterPerSecond, air_density: KilogramPerCubeMeter ) -> SquareMeterPerSecond: r"""Compute the kinematic viscosity of air. @@ -167,169 +164,6 @@ def compute_kinematic_viscosity_of_air( return dynamic_viscosity_of_air / air_density -# Unitless quantities -##################### - - -def compute_reynolds_number( - wind_speed: MeterPerSecond, - conductor_diameter: Meter, - kinematic_viscosity_of_air: SquareMeterPerSecond, -) -> Unitless: - r"""Compute the Reynolds number using the conductor diameter as characteristic length scale. - - Defined in the text on page 25 of :cite:p:`cigre601`. - - The Reynolds number is a dimensionless quantity that can be used to assess if a stream is - likely to be turbulent or not. It is given by - - .. math:: - - \text{Re} = \frac{v L}{\nu}, - - where :math:`v` is the flow velocity, :math:`L` is a *characteristic length* (in our case, - the conductor diameter) and :math:`\nu` is the kinematic viscosity. - - Parameters - ---------- - wind_speed: - :math:`v~\left[\text{m}~\text{s}^{-1}\right]`. The wind speed. - conductor_diameter: - :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. - kinematic_viscosity_of_air: - :math:`\nu_f~\left[\text{m}^2~\text{s}^{-1}\right]`. The kinematic viscosity of air. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`\text{Re}`. The Reynolds number. - """ - v = wind_speed - D = conductor_diameter - nu_f = kinematic_viscosity_of_air - return v * D / nu_f - - -def compute_grashof_number( - conductor_diameter: Meter, - surface_temperature: Celsius, - air_temperature: Celsius, - kinematic_viscosity_of_air: SquareMeterPerSecond, - coefficient_of_gravity: MeterPerSquareSecond = 9.807, -) -> Unitless: - r"""Compute the Grashof number. - - Defined in the nomenclature on page 7 of :cite:p:`cigre601`. - - The Grashof number is a dimensionless quantity that can be used to assess the degree of free - and forced convective heat transfer. - - Parameters - ---------- - conductor_diameter: - :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. - surface_temperature: - :math:`T_s~\left[^\circ\text{C}\right]`. The conductor surface temperature. - air_temperature: - :math:`T_a~\left[^\circ\text{C}\right]`. The ambient air temperature. - kinematic_viscosity_of_air: - :math:`\nu_f~\left[\text{m}^2~\text{s}^{-1}\right]`. The kinematic viscosity of air. - coefficient_of_gravity: - :math:`g~\left[\text{m}~\text{s}^{-2}\right]`. The graviatational constant, optional - (default=9.807). - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`\text{Gr}`. The Grashof number. - """ - T_s = surface_temperature - T_a = air_temperature - T_f = 0.5 * (T_s + T_a) - D = conductor_diameter - nu_f = kinematic_viscosity_of_air - g = coefficient_of_gravity - - return (D**3) * np.abs((T_s - T_a)) * g / ((T_f + 273.15) * (nu_f**2)) - - -def compute_prandtl_number( - thermal_conductivity_of_air: WattPerMeterPerKelvin, - dynamic_viscosity_of_air: KilogramPerMeterPerSecond, - specific_heat_capacity_of_air: JoulePerKilogramPerKelvin, -) -> Unitless: - r"""Compute the Prandtl number. - - Defined in the nomenclature on page 8 of :cite:p:`cigre601`. - - The Prandtl number measures the ratio between viscosity and thermal diffusivity for a fluid. - - Parameters - ---------- - thermal_conductivity_of_air: - :math:`\lambda_f~\left[\text{W}~\text{m}^{-1}~\text{K}^{-1}\right]`. The thermal - conductivity of air at the given temperature. - dynamic_viscosity_of_air: - :math:`\mu_f~\left[\text{kg}~\text{m}^{-1}~\text{s}^{-1}\right]`. The dynamic viscosity of - air. - specific_heat_capacity_of_air: - :math:`\text{J}~\left[\text{kg}^{-1}~\text{K}^{-1}\right]`. The specific heat capacity of - air. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`\text{Pr}`. The Prandtl number. - """ - lambda_f = thermal_conductivity_of_air - mu_f = dynamic_viscosity_of_air - c_f = specific_heat_capacity_of_air - - return c_f * mu_f / lambda_f - - -def compute_conductor_roughness( - conductor_diameter: Meter, - outer_layer_strand_diameter: Meter, -) -> Unitless: - r"""Compute the surface roughness of the conductor. - - Defined in the text on page 25 of :cite:p:`cigre601`. - - Parameters - ---------- - conductor_diameter: - :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. - outer_layer_strand_diameter: - :math:`d~\left[\text{m}\right]`. The diameter of the strands in the outer layer of the - conductor. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`\text{Rs}`. The roughness number - """ - D = conductor_diameter - d = outer_layer_strand_diameter - - if np.any(d < 0): - raise ValueError( - dedent( - """\ - Cannot have negative outer layer strand diameter. If it was set this way to signify - that it is a smooth conductor, then set the outer layer strand diameter to nan or zero - instead.\ - """ - ) - ) - if np.any(d >= D): - raise ValueError( - "The outer layer strand diameter must be strictly smaller than the conductor diameter." - ) - - return d / (2 * (D - d)) - - ## Nusselt number calculation ############################# @@ -351,8 +185,8 @@ def _check_perpendicular_flow_nusseltnumber_out_of_bounds(reynolds_number, condu @vectorize(nopython=True) def _compute_perpendicular_flow_nusseltnumber( - reynolds_number: Unitless, - conductor_roughness: Meter, + reynolds_number: Unitless, + conductor_roughness: Meter, ) -> Unitless: # TODO: Look at references for this table Re = reynolds_number @@ -382,12 +216,12 @@ def _compute_perpendicular_flow_nusseltnumber( else: B, n = 0.048, 0.800 - return B * Re**n # type: ignore + 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. @@ -419,9 +253,9 @@ def compute_perpendicular_flow_nusseltnumber( @vectorize(nopython=True) def _correct_wind_direction_effect_on_nusselt_number( - perpendicular_flow_nusselt_number: Unitless, - angle_of_attack: Radian, - conductor_roughness: Unitless, + perpendicular_flow_nusselt_number: Unitless, + angle_of_attack: Radian, + conductor_roughness: Unitless, ) -> Unitless: delta = angle_of_attack Nu_90 = perpendicular_flow_nusselt_number @@ -430,23 +264,23 @@ def _correct_wind_direction_effect_on_nusselt_number( sin_delta = np.sin(delta) if Rs == 0 or np.isnan(Rs): - sin_delta_sq = sin_delta**2 + sin_delta_sq = sin_delta ** 2 cos_delta_sq = 1 - sin_delta_sq correction_factor = (sin_delta_sq + cos_delta_sq * 0.0169) ** 0.225 else: if delta <= np.radians(24): - correction_factor = 0.42 + 0.68 * (sin_delta**1.08) + correction_factor = 0.42 + 0.68 * (sin_delta ** 1.08) else: - correction_factor = 0.42 + 0.58 * (sin_delta**0.90) + 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, - conductor_roughness: Unitless, + perpendicular_flow_nusselt_number: Unitless, + angle_of_attack: Radian, + conductor_roughness: Unitless, ) -> Unitless: r"""Correct the Nusselt number for the wind's angle-of-attack. @@ -481,7 +315,7 @@ def correct_wind_direction_effect_on_nusselt_number( 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): @@ -492,26 +326,26 @@ def _check_horizontal_natural_nusselt_number( @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 < 1e-1: return 0 elif GrPr < 1e2: - return 1.020 * GrPr**0.148 + return 1.020 * GrPr ** 0.148 elif GrPr < 1e4: - return 0.850 * GrPr**0.188 + return 0.850 * GrPr ** 0.188 elif GrPr < 1e7: - return 0.480 * GrPr**0.250 + return 0.480 * GrPr ** 0.250 else: - return 0.125 * GrPr**0.333 + return 0.125 * GrPr ** 0.333 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. @@ -544,8 +378,8 @@ def compute_horizontal_natural_nusselt_number( def _check_conductor_inclination( - conductor_inclination: Radian, - conductor_roughness: Unitless, + conductor_inclination: Radian, + conductor_roughness: Unitless, ) -> None: beta = np.degrees(conductor_inclination) Rs = conductor_roughness @@ -562,24 +396,24 @@ def _check_conductor_inclination( @vectorize(nopython=True) def _correct_natural_nusselt_number_inclination( - horizontal_natural_nusselt_number: Unitless, - conductor_inclination: Radian, - conductor_roughness: Unitless, + horizontal_natural_nusselt_number: Unitless, + conductor_inclination: Radian, + conductor_roughness: Unitless, ) -> Unitless: beta = np.degrees(conductor_inclination) Nu_nat = horizontal_natural_nusselt_number Rs = conductor_roughness if Rs == 0 or np.isnan(Rs): - return Nu_nat * (1 - 1.58e-4 * beta**1.5) + return Nu_nat * (1 - 1.58e-4 * beta ** 1.5) else: - return Nu_nat * (1 - 1.76e-6 * beta**2.5) + return Nu_nat * (1 - 1.76e-6 * beta ** 2.5) def correct_natural_nusselt_number_inclination( - horizontal_natural_nusselt_number: Unitless, - conductor_inclination: Radian, - conductor_roughness: Unitless, + horizontal_natural_nusselt_number: Unitless, + conductor_inclination: Radian, + conductor_roughness: Unitless, ) -> Unitless: r"""Correct the natural Nusselt number for the effect of the span inclination. @@ -613,8 +447,8 @@ def correct_natural_nusselt_number_inclination( def compute_nusselt_number( - forced_convection_nusselt_number: Unitless, - natural_nusselt_number: Unitless, + forced_convection_nusselt_number: Unitless, + natural_nusselt_number: Unitless, ) -> Unitless: r"""Compute the nusselt number. @@ -634,44 +468,3 @@ def compute_nusselt_number( :math:`Nu`. The nusselt number. """ return np.maximum(forced_convection_nusselt_number, natural_nusselt_number) - - -# Cooling computation -##################### - - -def compute_convective_cooling( - surface_temperature: Celsius, - air_temperature: Celsius, - nusselt_number: Unitless, - thermal_conductivity_of_air: WattPerMeterPerKelvin, -) -> WattPerMeter: - r"""Compute the convective cooling of the conductor. - - Equation (17) on page 24 of :cite:p:`cigre601`. - - Parameters - ---------- - surface_temperature: - :math:`T_s~\left[^\circ\text{C}\right]`. The conductor surface temperature. - air_temperature: - :math:`T_a~\left[^\circ\text{C}\right]`. The ambient air temperature. - nusselt_number: - :math:`Nu`. The nusselt number. - thermal_conductivity_of_air: - :math:`\lambda_f~\left[\text{W}~\text{m}^{-1}~\text{K}^{-1}\right]`. The thermal - conductivity of air at the given temperature. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_c~\left[\text{W}~\text{m}^{-1}\right]`. The convective cooling of the conductor. - Either due to wind, or passive convection, whichever is largest. - """ - pi = np.pi - lambda_f = thermal_conductivity_of_air - T_s = surface_temperature - T_a = air_temperature - Nu = nusselt_number - - return pi * lambda_f * (T_s - T_a) * Nu diff --git a/linerate/equations/convective_cooling.py b/linerate/equations/convective_cooling.py new file mode 100644 index 0000000..cec9a54 --- /dev/null +++ b/linerate/equations/convective_cooling.py @@ -0,0 +1,40 @@ +import numpy as np + +from linerate.units import Celsius, Unitless, WattPerMeterPerKelvin, WattPerMeter + + +def compute_convective_cooling( + surface_temperature: Celsius, + air_temperature: Celsius, + nusselt_number: Unitless, + thermal_conductivity_of_air: WattPerMeterPerKelvin, +) -> WattPerMeter: + r"""Compute the convective cooling of the conductor. + + Equation (17) on page 24 of :cite:p:`cigre601` and Equation (12) on page 6 of :cite:p:`cigre207`. + + Parameters + ---------- + surface_temperature: + :math:`T_s~\left[^\circ\text{C}\right]`. The conductor surface temperature. + air_temperature: + :math:`T_a~\left[^\circ\text{C}\right]`. The ambient air temperature. + nusselt_number: + :math:`Nu`. The nusselt number. + thermal_conductivity_of_air: + :math:`\lambda_f~\left[\text{W}~\text{m}^{-1}~\text{K}^{-1}\right]`. The thermal + conductivity of air at the given temperature. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_c~\left[\text{W}~\text{m}^{-1}\right]`. The convective cooling of the conductor. + Either due to wind, or passive convection, whichever is largest. + """ + pi = np.pi + lambda_f = thermal_conductivity_of_air + T_s = surface_temperature + T_a = air_temperature + Nu = nusselt_number + + return pi * lambda_f * (T_s - T_a) * Nu \ No newline at end of file diff --git a/linerate/equations/dimensionless.py b/linerate/equations/dimensionless.py new file mode 100644 index 0000000..10430ad --- /dev/null +++ b/linerate/equations/dimensionless.py @@ -0,0 +1,166 @@ +from textwrap import dedent + +import numpy as np + +from linerate.units import Unitless, WattPerMeterPerKelvin, KilogramPerMeterPerSecond, JoulePerKilogramPerKelvin, \ + MeterPerSecond, Meter, SquareMeterPerSecond, Celsius, MeterPerSquareSecond + + +def compute_reynolds_number( + wind_speed: MeterPerSecond, + conductor_diameter: Meter, + kinematic_viscosity_of_air: SquareMeterPerSecond, +) -> Unitless: + r"""Compute the Reynolds number using the conductor diameter as characteristic length scale. + + Defined in the text on page 25 of :cite:p:`cigre601` and equation (2c) on page 10 in :cite:p:`ieee738`. + + The Reynolds number is a dimensionless quantity that can be used to assess if a stream is + likely to be turbulent or not. It is given by + + .. math:: + + \text{Re} = \frac{v L}{\nu}, + + where :math:`v` is the flow velocity, :math:`L` is a *characteristic length* (in our case, + the conductor diameter) and :math:`\nu` is the kinematic viscosity. + + Parameters + ---------- + wind_speed: + :math:`v~\left[\text{m}~\text{s}^{-1}\right]`. The wind speed. + conductor_diameter: + :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. + kinematic_viscosity_of_air: + :math:`\nu_f~\left[\text{m}^2~\text{s}^{-1}\right]`. The kinematic viscosity of air. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`\text{Re}`. The Reynolds number. + """ + v = wind_speed + D = conductor_diameter + nu_f = kinematic_viscosity_of_air + return v * D / nu_f + + +def compute_grashof_number( + conductor_diameter: Meter, + surface_temperature: Celsius, + air_temperature: Celsius, + kinematic_viscosity_of_air: SquareMeterPerSecond, + coefficient_of_gravity: MeterPerSquareSecond = 9.807, +) -> Unitless: + r"""Compute the Grashof number. + + Defined in the nomenclature on page 7 of :cite:p:`cigre601` + and on page 5 of :cite:p:`cigre207`. + + The Grashof number is a dimensionless quantity that can be used to assess the degree of free + and forced convective heat transfer. + + Parameters + ---------- + conductor_diameter: + :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. + surface_temperature: + :math:`T_s~\left[^\circ\text{C}\right]`. The conductor surface temperature. + air_temperature: + :math:`T_a~\left[^\circ\text{C}\right]`. The ambient air temperature. + kinematic_viscosity_of_air: + :math:`\nu_f~\left[\text{m}^2~\text{s}^{-1}\right]`. The kinematic viscosity of air. + coefficient_of_gravity: + :math:`g~\left[\text{m}~\text{s}^{-2}\right]`. The graviatational constant, optional + (default=9.807). + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`\text{Gr}`. The Grashof number. + """ + T_s = surface_temperature + T_a = air_temperature + T_f = 0.5 * (T_s + T_a) + D = conductor_diameter + nu_f = kinematic_viscosity_of_air + g = coefficient_of_gravity + + return (D**3) * np.abs((T_s - T_a)) * g / ((T_f + 273.15) * (nu_f**2)) + + +def compute_prandtl_number( + thermal_conductivity_of_air: WattPerMeterPerKelvin, + dynamic_viscosity_of_air: KilogramPerMeterPerSecond, + specific_heat_capacity_of_air: JoulePerKilogramPerKelvin, +) -> Unitless: + r"""Compute the Prandtl number. + + Defined in the nomenclature on page 8 of :cite:p:`cigre601`. + + The Prandtl number measures the ratio between viscosity and thermal diffusivity for a fluid. + + Parameters + ---------- + thermal_conductivity_of_air: + :math:`\lambda_f~\left[\text{W}~\text{m}^{-1}~\text{K}^{-1}\right]`. The thermal + conductivity of air at the given temperature. + dynamic_viscosity_of_air: + :math:`\mu_f~\left[\text{kg}~\text{m}^{-1}~\text{s}^{-1}\right]`. The dynamic viscosity of + air. + specific_heat_capacity_of_air: + :math:`\text{J}~\left[\text{kg}^{-1}~\text{K}^{-1}\right]`. The specific heat capacity of + air. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`\text{Pr}`. The Prandtl number. + """ + lambda_f = thermal_conductivity_of_air + mu_f = dynamic_viscosity_of_air + c_f = specific_heat_capacity_of_air + + return c_f * mu_f / lambda_f + + +def compute_conductor_roughness( + conductor_diameter: Meter, + outer_layer_strand_diameter: Meter, +) -> Unitless: + r"""Compute the surface roughness of the conductor. + + Defined in the text on page 25 of :cite:p:`cigre601` and on page 6 of :cite:p:`cigre207`. + + Parameters + ---------- + conductor_diameter: + :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. + outer_layer_strand_diameter: + :math:`d~\left[\text{m}\right]`. The diameter of the strands in the outer layer of the + conductor. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`\text{Rs}`. The roughness number + """ + D = conductor_diameter + d = outer_layer_strand_diameter + + if np.any(d < 0): + raise ValueError( + dedent( + """\ + Cannot have negative outer layer strand diameter. If it was set this way to signify + that it is a smooth conductor, then set the outer layer strand diameter to nan or zero + instead.\ + """ + ) + ) + if np.any(d >= D): + raise ValueError( + "The outer layer strand diameter must be strictly smaller than the conductor diameter." + ) + + return d / (2 * (D - d)) diff --git a/linerate/equations/ieee738/convective_cooling.py b/linerate/equations/ieee738/convective_cooling.py index 33cf62c..aac5e75 100644 --- a/linerate/equations/ieee738/convective_cooling.py +++ b/linerate/equations/ieee738/convective_cooling.py @@ -5,7 +5,6 @@ KilogramPerCubeMeter, KilogramPerMeterPerSecond, Meter, - MeterPerSecond, Radian, SquareMeterPerSecond, Unitless, @@ -91,50 +90,6 @@ def compute_kinematic_viscosity_of_air( # nu_f return dynamic_viscosity_of_air / air_density -def compute_reynolds_number( - wind_speed: MeterPerSecond, - conductor_diameter: Meter, - kinematic_viscosity_of_air: SquareMeterPerSecond, -) -> Unitless: - r"""Compute the Reynolds number using the conductor diameter as characteristic length scale. - - Equation (2c) on page 10 in :cite:p:`ieee738` - - Separated out kinematic viscosity of air into its own function because that is how CIGRE601 - has done it, and that way they do the calculation the same way. Kinematic viscosity is given by - - .. math: - - \nu_f = \frac{\mu_f}{\rho_f} - - where {\mu_f} is the dynamic viscosity of air and {\rho_f} is the air density. - - The Reynolds number is a dimensionless quantity that can be used to assess if a stream is - likely to be turbulent or not. It is given by - - .. math:: - - \text{N_{Re}} = \frac{V_w D_0}{\nu_f}, - - where :math:`V_w` is the flow velocity, :math:`D_0` is a *characteristic length* (in our case, - the conductor diameter) and :math:`\nu_f` is the kinematic viscosity. - - Parameters - ---------- - wind_speed: - :math:`V_w~\left[\text{m}~\text{s}^{-1}\right]`. The wind speed. - conductor_diameter: - :math:`D_0~\left[\text{m}\right]`. Outer diameter of the conductor. - kinematic_viscosity_of_air: - :math:`\nu_f~\left[\text{m}^2~\text{s}^{-1}\right]`. The kinematic viscosity of air. - - """ - V_w = wind_speed - D_0 = conductor_diameter - nu_f = kinematic_viscosity_of_air - return V_w * D_0 / nu_f - - def compute_wind_direction_factor( # K_angle angle_of_attack: Radian, ) -> Unitless: diff --git a/linerate/model.py b/linerate/model.py index 3977c18..3543298 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -19,7 +19,7 @@ joule_heating, math, radiative_cooling, - solar_angles, + solar_angles, convective_cooling, dimensionless, ) from linerate.equations.math import switch_cos_sin from linerate.types import Span, Weather @@ -386,12 +386,12 @@ def compute_convective_cooling( # Compute unitless quantities Re = np.minimum( - cigre601.convective_cooling.compute_reynolds_number(V, D, nu_f), + dimensionless.compute_reynolds_number(V, D, nu_f), self.max_reynolds_number, ) - Gr = cigre601.convective_cooling.compute_grashof_number(D, T_c, T_a, nu_f) - Pr = cigre601.convective_cooling.compute_prandtl_number(lambda_f, mu_f, c_f) - Rs = cigre601.convective_cooling.compute_conductor_roughness(D, d) + Gr = dimensionless.compute_grashof_number(D, T_c, T_a, nu_f) + Pr = dimensionless.compute_prandtl_number(lambda_f, mu_f, c_f) + Rs = dimensionless.compute_conductor_roughness(D, d) # Compute nusselt numbers Nu_90 = cigre601.convective_cooling.compute_perpendicular_flow_nusseltnumber( @@ -410,7 +410,7 @@ def compute_convective_cooling( forced_convection_nusselt_number=Nu_delta, natural_nusselt_number=Nu_beta ) - return cigre601.convective_cooling.compute_convective_cooling( + return convective_cooling.compute_convective_cooling( surface_temperature=conductor_temperature, air_temperature=self.weather.air_temperature, nusselt_number=Nu, @@ -519,7 +519,7 @@ def compute_convective_cooling( rho_f = ieee738.convective_cooling.compute_air_density(T_f, y) nu_f = ieee738.convective_cooling.compute_kinematic_viscosity_of_air(mu_f, rho_f) Re = np.minimum( - ieee738.convective_cooling.compute_reynolds_number(V, D, nu_f), # N_Re in IEEE + dimensionless.compute_reynolds_number(V, D, nu_f), # N_Re in IEEE self.max_reynolds_number, ) delta = math.compute_angle_of_attack( diff --git a/linerate/models/Cigre207.md b/linerate/models/Cigre207.md new file mode 100644 index 0000000..b3702ed --- /dev/null +++ b/linerate/models/Cigre207.md @@ -0,0 +1,30 @@ +# Cigre 207 + +## Cigre 207 vs 601 + +The Cigre 601 report states that 601 takes into account improvements in the calculation of AC resistance of stranded +conductors, includes a more realistic estimation of the axial and radial temeprature distributions, more coherent +convection model for low wind speeds, and has a more flexible solar radiation model. + +## Skin effect + +$P_J = k_j I^2 R_{dc}$ +Where $k_j$ is the skin effect factor. +We do not implement skin effect, since this is included in the AC resistance values we use as input. + +## Magnetic effect + +The effect of magnetic heating seems to be intertwined with the skin effect in Cigre 207. +To make things simpler and hopefully accurate enough, we use the same linearised model as our implementation of +Cigre 601. +**Maybe we have to come back to this one...** + +## Solar heating + +Solar heating is the same as for Cigre 601, but with clearness ratio 1. + +## Convective cooling +There are a few inconsistencies in the Cigre 207 modelling of convective cooling. +The definition of the Reynolds number is given as +$Re = \rho_r V D/\nu_f$, where \rho_r is the air density relative to density at sea level. +This factor is not present in the ordinary definition of Reynolds number, so we choose to omit it. \ No newline at end of file diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py new file mode 100644 index 0000000..c15c0e3 --- /dev/null +++ b/linerate/models/cigre207.py @@ -0,0 +1,117 @@ +from linerate import ThermalModel, Span, Weather +from linerate.equations import solar_angles, cigre601, math, cigre207, dimensionless, convective_cooling +from linerate.equations.math import switch_cos_sin +from linerate.model import _copy_method_docstring +from linerate.units import Date, Celsius, Ampere, WattPerMeter + + +class Cigre207(ThermalModel): + def __init__( + self, + span: Span, + weather: Weather, + time: Date, + ): + super().__init__(span, weather) + self.time = time + + @_copy_method_docstring(ThermalModel) + def compute_joule_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + return super().compute_joule_heating( + conductor_temperature=conductor_temperature, current=current + ) + + @_copy_method_docstring(ThermalModel) + def compute_solar_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + alpha_s = self.span.conductor.solar_absorptivity + F = self.span.ground_albedo + phi = self.span.latitude + gamma_c = self.span.conductor_azimuth + y = self.span.conductor_altitude + N_s = self.weather.clearness_ratio + D = self.span.conductor.conductor_diameter + + omega = solar_angles.compute_hour_angle_relative_to_noon(self.time) + delta = solar_angles.compute_solar_declination(self.time) + sin_H_s = solar_angles.compute_sin_solar_altitude(phi, delta, omega) + chi = solar_angles.compute_solar_azimuth_variable(phi, delta, omega) + C = solar_angles.compute_solar_azimuth_constant(chi, omega) + gamma_s = solar_angles.compute_solar_azimuth(C, chi) # Z_c in IEEE + cos_eta = solar_angles.compute_cos_solar_effective_incidence_angle( + sin_H_s, gamma_s, gamma_c + ) + sin_eta = switch_cos_sin(cos_eta) + + I_B = cigre601.solar_heating.compute_direct_solar_radiation(sin_H_s, N_s, y) + I_d = cigre601.solar_heating.compute_diffuse_sky_radiation(I_B, sin_H_s) + I_T = cigre601.solar_heating.compute_global_radiation_intensity( + I_B, I_d, F, sin_eta, sin_H_s + ) + return cigre601.solar_heating.compute_solar_heating( + alpha_s, + I_T, + D, + ) + + @_copy_method_docstring(ThermalModel) + def compute_convective_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + D = self.span.conductor.conductor_diameter + d = self.span.conductor.outer_layer_strand_diameter + y = self.span.conductor_altitude + V = self.weather.wind_speed + T_a = self.weather.air_temperature + T_c = conductor_temperature + T_f = 0.5 * (T_c + T_a) + + # Compute physical quantities + lambda_f = cigre207.convective_cooling.compute_thermal_conductivity_of_air(T_f) + nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) + delta = math.compute_angle_of_attack( + self.weather.wind_direction, self.span.conductor_azimuth + ) + + # Compute unitless quantities + # Reynolds number is defined in the text on page 5 of :cite:p:`cigre207`. + # The definition includes a relative air density, which does not make sense, so we omit it here and use the + # standard definition of Reynolds number instead. + Re = dimensionless.compute_reynolds_number(V, D, nu_f) + Gr = dimensionless.compute_grashof_number(D, T_c, T_a, nu_f) + Pr = cigre207.convective_cooling.compute_prandtl_number(T_f) + Rs = dimensionless.compute_conductor_roughness(D, d) + + # Compute nusselt numbers + Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber( + reynolds_number=Re, conductor_roughness=Rs + ) + Nu_delta = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number( + Nu_90, delta + ) + Nu_cor = cigre207.convective_cooling.compute_low_wind_speed_nusseltnumber(Nu_90) + + Nu_0 = cigre207.convective_cooling.compute_horizontal_natural_nusselt_number(Gr, Pr) + + Nu = cigre207.convective_cooling.compute_nusselt_number( + forced_convection_nusselt_number=Nu_delta, natural_nusselt_number=Nu_0, + low_wind_nusselt_number=Nu_cor, wind_speed=V + ) + + return convective_cooling.compute_convective_cooling( + surface_temperature=conductor_temperature, + air_temperature=self.weather.air_temperature, + nusselt_number=Nu, + thermal_conductivity_of_air=lambda_f, + ) + + @_copy_method_docstring(ThermalModel) + def compute_radiative_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + return super().compute_radiative_cooling( + conductor_temperature=conductor_temperature, current=current + ) From 7202959dbb8b7e4855433813c37d122e4fca23e9 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Wed, 29 May 2024 14:20:09 +0200 Subject: [PATCH 02/18] Fix tests, compute Reynolds number as in Cigre207 report --- linerate/equations/cigre207/ac_resistance.py | 34 +-- .../equations/cigre207/convective_cooling.py | 38 ++- linerate/models/Cigre207.md | 2 +- linerate/models/cigre207.py | 49 +++- .../cigre601/test_convective_cooling.py | 241 +----------------- ...ing.py => test_ieee_convective_cooling.py} | 36 --- tests/equations/test_convective_cooling.py | 97 +++++++ tests/equations/test_dimensionless.py | 157 ++++++++++++ 8 files changed, 346 insertions(+), 308 deletions(-) rename tests/equations/ieee738/{test_convective_cooling.py => test_ieee_convective_cooling.py} (89%) create mode 100644 tests/equations/test_convective_cooling.py create mode 100644 tests/equations/test_dimensionless.py diff --git a/linerate/equations/cigre207/ac_resistance.py b/linerate/equations/cigre207/ac_resistance.py index 86fa2ff..748baef 100644 --- a/linerate/equations/cigre207/ac_resistance.py +++ b/linerate/equations/cigre207/ac_resistance.py @@ -1,45 +1,27 @@ from typing import Union -from linerate.units import OhmPerMeter, Ampere, SquareMeter, Unitless, SquareMeterPerAmpere +from linerate.units import OhmPerMeter -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, +def correct_resistance_for_skin_effect( + dc_resistance: OhmPerMeter, ) -> OhmPerMeter: r""" - Return resistance with constant correction for magnetic effects, using simple method from + Return resistance with constant correction for skin effect, 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. + dc_resistance: + :math:`R~\left[\Omega\right]`. The DC resistance of the conductor. 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. + taking skin effect into account. """ - return 1.0123 * ac_resistance + return 1.0123 * dc_resistance # TODO: Implement section 2.1.2? \ No newline at end of file diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py index e3b12f7..78e7215 100644 --- a/linerate/equations/cigre207/convective_cooling.py +++ b/linerate/equations/cigre207/convective_cooling.py @@ -10,7 +10,7 @@ MeterPerSecond, Radian, Unitless, - WattPerMeterPerKelvin, + WattPerMeterPerKelvin, SquareMeterPerSecond, ) @@ -107,6 +107,42 @@ def compute_prandtl_number( return 0.715 - 2.5e-4*film_temperature +def compute_reynolds_number( + 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. + + Defined on page 5 of :cite:p:`cigre207`. + This is a non-standard definition which seems to indicate that the kinematic viscosity has to + be corrected for the density. + + Parameters + ---------- + wind_speed: + :math:`v~\left[\text{m}~\text{s}^{-1}\right]`. The wind speed. + conductor_diameter: + :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. + kinematic_viscosity_of_air: + :math:`\nu_f~\left[\text{m}^2~\text{s}^{-1}\right]`. The kinematic viscosity of air. + relative_air_density: + :math:`\rho_r~1`. The air density relative to density at sea level. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`\text{Re}`. The Reynolds number. + """ + v = wind_speed + D = conductor_diameter + nu_f = kinematic_viscosity_of_air + rho_r = relative_air_density + return rho_r * v * D / nu_f + + + ## Nusselt number calculation ############################# diff --git a/linerate/models/Cigre207.md b/linerate/models/Cigre207.md index b3702ed..3912193 100644 --- a/linerate/models/Cigre207.md +++ b/linerate/models/Cigre207.md @@ -3,7 +3,7 @@ ## Cigre 207 vs 601 The Cigre 601 report states that 601 takes into account improvements in the calculation of AC resistance of stranded -conductors, includes a more realistic estimation of the axial and radial temeprature distributions, more coherent +conductors, includes a more realistic estimation of the axial and radial temperature distributions, more coherent convection model for low wind speeds, and has a more flexible solar radiation model. ## Skin effect diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index c15c0e3..9f2ae92 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -1,8 +1,10 @@ +from abc import abstractmethod + from linerate import ThermalModel, Span, Weather -from linerate.equations import solar_angles, cigre601, math, cigre207, dimensionless, convective_cooling +from linerate.equations import solar_angles, cigre601, math, cigre207, dimensionless, convective_cooling, joule_heating from linerate.equations.math import switch_cos_sin from linerate.model import _copy_method_docstring -from linerate.units import Date, Celsius, Ampere, WattPerMeter +from linerate.units import Date, Celsius, Ampere, WattPerMeter, OhmPerMeter class Cigre207(ThermalModel): @@ -63,8 +65,8 @@ def compute_convective_cooling( ) -> WattPerMeter: D = self.span.conductor.conductor_diameter d = self.span.conductor.outer_layer_strand_diameter - y = self.span.conductor_altitude V = self.weather.wind_speed + y = self.span.conductor_altitude T_a = self.weather.air_temperature T_c = conductor_temperature T_f = 0.5 * (T_c + T_a) @@ -80,7 +82,8 @@ def compute_convective_cooling( # Reynolds number is defined in the text on page 5 of :cite:p:`cigre207`. # The definition includes a relative air density, which does not make sense, so we omit it here and use the # standard definition of Reynolds number instead. - Re = dimensionless.compute_reynolds_number(V, D, nu_f) + rho_r = cigre207.convective_cooling.compute_relative_air_density(y) + Re = cigre207.convective_cooling.compute_reynolds_number(V, D, nu_f, rho_r) Gr = dimensionless.compute_grashof_number(D, T_c, T_a, nu_f) Pr = cigre207.convective_cooling.compute_prandtl_number(T_f) Rs = dimensionless.compute_conductor_roughness(D, d) @@ -115,3 +118,41 @@ def compute_radiative_cooling( return super().compute_radiative_cooling( conductor_temperature=conductor_temperature, current=current ) + + @abstractmethod + def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: + r"""Compute the conductor resistance, :math:`R~\left[\Omega~\text{m}^{-1}\right]`. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`R~\left[\Omega\right]`. The resistance at the given temperature and current. + """ + resistance = joule_heating.compute_resistance( + conductor_temperature, + temperature1=self.span.conductor.temperature1, + temperature2=self.span.conductor.temperature2, + resistance_at_temperature1=self.span.conductor.resistance_at_temperature1, + resistance_at_temperature2=self.span.conductor.resistance_at_temperature2, + ) + + A = self.span.conductor.aluminium_cross_section_area + b = self.span.conductor.constant_magnetic_effect + m = self.span.conductor.current_density_proportional_magnetic_effect + max_increase = self.span.conductor.max_magnetic_core_relative_resistance_increase + + return joule_heating.correct_resistance_acsr_magnetic_core_loss( + ac_resistance=resistance, + current=current, + aluminium_cross_section_area=A, + constant_magnetic_effect=b, + current_density_proportional_magnetic_effect=m, + max_relative_increase=max_increase, + ) diff --git a/tests/equations/cigre601/test_convective_cooling.py b/tests/equations/cigre601/test_convective_cooling.py index 13a74c4..6e690c4 100644 --- a/tests/equations/cigre601/test_convective_cooling.py +++ b/tests/equations/cigre601/test_convective_cooling.py @@ -7,6 +7,7 @@ import linerate.equations.cigre601.convective_cooling as convective_cooling + # Tests for physical quantities ############################### @@ -259,155 +260,6 @@ def test_kinematic_viscosity_of_air_with_example(): assert convective_cooling.compute_kinematic_viscosity_of_air(mu_f, gamma) == approx(1) -# Tests for unitless quantities -############################### - - -@hypothesis.given(wind_speed=st.floats(min_value=0, max_value=1000, allow_nan=False)) -def test_reynolds_number_scales_linearly_with_wind_speed(wind_speed): - V = wind_speed - D = 1 - Nu_f = 1 - - assert convective_cooling.compute_reynolds_number(V, D, Nu_f) == pytest.approx(V) - - -@hypothesis.given(conductor_diameter=st.floats(min_value=0, max_value=1000, allow_nan=False)) -def test_reynolds_number_scales_linearly_with_diameter(conductor_diameter): - V = 1 - D = conductor_diameter - Nu_f = 1 - - assert convective_cooling.compute_reynolds_number(V, D, Nu_f) == pytest.approx(D) - - -@hypothesis.given(kinematic_viscosity=st.floats(min_value=1e-10, allow_nan=False)) -def test_reynolds_number_scales_inversly_with_kinematic_viscosity(kinematic_viscosity): - V = 1 - D = 1 - Nu_f = kinematic_viscosity - - assert convective_cooling.compute_reynolds_number(V, D, Nu_f) == approx(1 / Nu_f) - - -def test_reynolds_number_with_example(): - V = 0.1 - D = 1.2 - Nu_f = 10 - - assert convective_cooling.compute_reynolds_number(V, D, Nu_f) == approx(0.012) - - -@hypothesis.given(conductor_diameter=st.floats(min_value=1e-5, max_value=1e5, allow_nan=False)) -def test_grashof_number_scales_cubic_with_conductor_diameter(conductor_diameter): - D = conductor_diameter - T_s = 1 - T_a = 0 - Nu_f = 1 - g = 273.65 # 0.5*T_s in Kelvin - - Gr = convective_cooling.compute_grashof_number(D, T_s, T_a, Nu_f, g) - assert Gr == approx(D**3) - - -@hypothesis.given(surface_temperature=st.floats(min_value=273, max_value=1000, allow_nan=False)) -def test_grashof_number_scales_correctly_with_surface_temperature(surface_temperature): - D = 1 - T_s = surface_temperature - T_a = 0 - Nu_f = 1 - g = 0.5 * T_s + 273.15 - - Gr = convective_cooling.compute_grashof_number(D, T_s, T_a, Nu_f, g) - assert Gr == approx(abs(T_s)) - - -@hypothesis.given(air_temperature=st.floats(min_value=273, max_value=1000, allow_nan=False)) -def test_grashof_number_scales_correctly_with_air_temperature(air_temperature): - D = 1 - T_s = 0 - T_a = air_temperature - Nu_f = 1 - g = 0.5 * T_a + 273.15 - - Gr = convective_cooling.compute_grashof_number(D, T_s, T_a, Nu_f, g) - assert Gr == approx(abs(T_a)) - - -@hypothesis.given( - kinematic_viscosity_of_air=st.floats(min_value=1e-5, max_value=1e10, allow_nan=False) -) -def test_grashof_number_scales_inversely_squared_with_kinematic_viscosity( - kinematic_viscosity_of_air, -): - D = 1 - T_s = 1 - T_a = 0 - Nu_f = kinematic_viscosity_of_air - g = 273.65 # 0.5*T_s in Kelvin - - Gr = convective_cooling.compute_grashof_number(D, T_s, T_a, Nu_f, g) - assert Gr == approx(1 / (Nu_f**2)) - - -@hypothesis.given(coefficient_of_gravity=st.floats(min_value=1e-3, max_value=1e3, allow_nan=False)) -def test_grashof_number_scales_linearly_with_gravity(coefficient_of_gravity): - D = 1 - T_s = 1 - T_a = 0 - Nu_f = np.sqrt(1 / 273.65) - g = coefficient_of_gravity - - Gr = convective_cooling.compute_grashof_number(D, T_s, T_a, Nu_f, g) - assert Gr == approx(g) - - -def test_grashof_number_with_example(): - D = 4 - T_s = 4 - T_a = 2 - Nu_f = 2 - g = 276.15 - - Gr = convective_cooling.compute_grashof_number(D, T_s, T_a, Nu_f, g) - assert Gr == approx(32) - - -@hypothesis.given(thermal_conductivity=st.floats(min_value=1e-10, max_value=1e10, allow_nan=False)) -def test_prandtl_number_scales_inversely_with_thermal_conductivity(thermal_conductivity): - c_f = 1 - mu_f = 1 - lambda_f = thermal_conductivity - - assert convective_cooling.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(1 / lambda_f) - - -@hypothesis.given(dynamic_viscosity=st.floats(min_value=1e-10, max_value=1e10, allow_nan=False)) -def test_prandtl_number_scales_linearly_with_dynamic_viscosity(dynamic_viscosity): - c_f = 1 - mu_f = dynamic_viscosity - lambda_f = 1 - - assert convective_cooling.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(mu_f) - - -@hypothesis.given( - specific_heat_capacity=st.floats(min_value=1e-10, max_value=1e10, allow_nan=False) -) -def test_prandtl_number_scales_linearly_with_specific_heat_capacity(specific_heat_capacity): - c_f = specific_heat_capacity - mu_f = 1 - lambda_f = 1 - - assert convective_cooling.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(c_f) - - -def test_compute_prandtl_number_with_example(): - c_f = 0.5 - mu_f = 3 - lambda_f = 0.5 - - assert convective_cooling.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(3) ## Nusselt number calculation @@ -646,94 +498,3 @@ def test_compute_nusselt_number(forced_convection_nusselt_number, natural_nussel Nu_est = convective_cooling.compute_nusselt_number(Nu_delta, Nu_beta) assert Nu == Nu_est - - -@hypothesis.given( - nusselt_number=st.floats( - min_value=0, - max_value=1e10, - allow_nan=False, - ) -) -def test_convective_cooling_scales_linearly_with_nusselt_number(nusselt_number): - Nu = nusselt_number - T_s = 1 - T_a = 0 - lambda_f = 1 / np.pi - - assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx(Nu) - - -@hypothesis.given( - thermal_conductivity_of_air=st.floats( - min_value=0, - max_value=1e10, - allow_nan=False, - ) -) -def test_convective_cooling_scales_linearly_with_thermal_conductivity_of_air( - thermal_conductivity_of_air, -): - Nu = 1 / np.pi - T_s = 1 - T_a = 0 - lambda_f = thermal_conductivity_of_air - - assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx(lambda_f) - - -@hypothesis.given( - surface_temperature=st.floats( - min_value=-1e10, - max_value=1e10, - allow_nan=False, - ), - air_temperature=st.floats( - min_value=-1e10, - max_value=1e10, - allow_nan=False, - ), -) -def test_convective_cooling_scales_linearly_with_temperature_difference( - surface_temperature, air_temperature -): - Nu = 1 - T_s = surface_temperature - T_a = air_temperature - lambda_f = 1 / np.pi - - assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx( - T_s - T_a - ) - - -@hypothesis.given( - air_temperature=st.floats( - min_value=-1e10, - max_value=1e10, - allow_nan=False, - ) -) -def test_convective_cooling_scales_affinely_with_air_temperature(air_temperature): - Nu = 1 - T_s = 1 - T_a = air_temperature - lambda_f = 1 / np.pi - - assert 1 - convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx(T_a) - - -@hypothesis.given( - surface_temperature=st.floats( - min_value=-1e10, - max_value=1e10, - allow_nan=False, - ), -) -def test_convective_cooling_scales_affinely_with_surface_temperature(surface_temperature): - Nu = 1 - T_s = surface_temperature - T_a = 1 - lambda_f = 1 / np.pi - - assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) + 1 == approx(T_s) diff --git a/tests/equations/ieee738/test_convective_cooling.py b/tests/equations/ieee738/test_ieee_convective_cooling.py similarity index 89% rename from tests/equations/ieee738/test_convective_cooling.py rename to tests/equations/ieee738/test_ieee_convective_cooling.py index 2d9a887..2d7c97e 100644 --- a/tests/equations/ieee738/test_convective_cooling.py +++ b/tests/equations/ieee738/test_ieee_convective_cooling.py @@ -1,7 +1,6 @@ import hypothesis import hypothesis.strategies as st import numpy as np -import pytest from pytest import approx from scipy.interpolate import lagrange @@ -102,41 +101,6 @@ def test_kinematic_viscosity_of_air_with_example(): assert convective_cooling.compute_kinematic_viscosity_of_air(mu_f, rho_f) == approx(1) -@hypothesis.given(wind_speed=st.floats(min_value=0, max_value=1000, allow_nan=False)) -def test_reynolds_number_scales_linearly_with_wind_speed(wind_speed): - v = wind_speed - D = 1 - nu_f = 1 - - assert convective_cooling.compute_reynolds_number(v, D, nu_f) == pytest.approx(v) - - -@hypothesis.given(conductor_diameter=st.floats(min_value=0, max_value=1000, allow_nan=False)) -def test_reynolds_number_scales_linearly_with_diameter(conductor_diameter): - v = 1 - D = conductor_diameter - nu_f = 1 - - assert convective_cooling.compute_reynolds_number(v, D, nu_f) == pytest.approx(D) - - -@hypothesis.given(kinematic_viscosity=st.floats(min_value=1e-10, allow_nan=False)) -def test_reynolds_number_scales_inversly_with_kinematic_viscosity(kinematic_viscosity): - v = 1 - D = 1 - nu_f = kinematic_viscosity - - assert convective_cooling.compute_reynolds_number(v, D, nu_f) == approx(1 / nu_f) - - -def test_reynolds_number_with_example(): - v = 0.1 - D = 1.2 - nu_f = 10 - - assert convective_cooling.compute_reynolds_number(v, D, nu_f) == approx(0.012) - - def test_wind_direction_factor_with_example(): Phi = 0 assert convective_cooling.compute_wind_direction_factor(Phi) == approx(0.388) diff --git a/tests/equations/test_convective_cooling.py b/tests/equations/test_convective_cooling.py new file mode 100644 index 0000000..876f749 --- /dev/null +++ b/tests/equations/test_convective_cooling.py @@ -0,0 +1,97 @@ +import hypothesis +import hypothesis.strategies as st +import numpy as np + +from pytest import approx +from linerate.equations import convective_cooling + + +@hypothesis.given( + nusselt_number=st.floats( + min_value=0, + max_value=1e10, + allow_nan=False, + ) +) +def test_convective_cooling_scales_linearly_with_nusselt_number(nusselt_number): + Nu = nusselt_number + T_s = 1 + T_a = 0 + lambda_f = 1 / np.pi + + assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx(Nu) + + +@hypothesis.given( + thermal_conductivity_of_air=st.floats( + min_value=0, + max_value=1e10, + allow_nan=False, + ) +) +def test_convective_cooling_scales_linearly_with_thermal_conductivity_of_air( + thermal_conductivity_of_air, +): + Nu = 1 / np.pi + T_s = 1 + T_a = 0 + lambda_f = thermal_conductivity_of_air + + assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx(lambda_f) + + +@hypothesis.given( + surface_temperature=st.floats( + min_value=-1e10, + max_value=1e10, + allow_nan=False, + ), + air_temperature=st.floats( + min_value=-1e10, + max_value=1e10, + allow_nan=False, + ), +) +def test_convective_cooling_scales_linearly_with_temperature_difference( + surface_temperature, air_temperature +): + Nu = 1 + T_s = surface_temperature + T_a = air_temperature + lambda_f = 1 / np.pi + + assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx( + T_s - T_a + ) + + +@hypothesis.given( + air_temperature=st.floats( + min_value=-1e10, + max_value=1e10, + allow_nan=False, + ) +) +def test_convective_cooling_scales_affinely_with_air_temperature(air_temperature): + Nu = 1 + T_s = 1 + T_a = air_temperature + lambda_f = 1 / np.pi + + assert 1 - convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) == approx(T_a) + + +@hypothesis.given( + surface_temperature=st.floats( + min_value=-1e10, + max_value=1e10, + allow_nan=False, + ), +) +def test_convective_cooling_scales_affinely_with_surface_temperature(surface_temperature): + Nu = 1 + T_s = surface_temperature + T_a = 1 + lambda_f = 1 / np.pi + + assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) + 1 == approx(T_s) diff --git a/tests/equations/test_dimensionless.py b/tests/equations/test_dimensionless.py new file mode 100644 index 0000000..efbf4ba --- /dev/null +++ b/tests/equations/test_dimensionless.py @@ -0,0 +1,157 @@ +# Tests for unitless quantities +############################### + +import hypothesis +import hypothesis.strategies as st +import numpy as np +import pytest +from pytest import approx + +from linerate.equations import dimensionless + + +@hypothesis.given(wind_speed=st.floats(min_value=0, max_value=1000, allow_nan=False)) +def test_reynolds_number_scales_linearly_with_wind_speed(wind_speed): + V = wind_speed + D = 1 + Nu_f = 1 + + assert dimensionless.compute_reynolds_number(V, D, Nu_f) == pytest.approx(V) + + +@hypothesis.given(conductor_diameter=st.floats(min_value=0, max_value=1000, allow_nan=False)) +def test_reynolds_number_scales_linearly_with_diameter(conductor_diameter): + V = 1 + D = conductor_diameter + Nu_f = 1 + + assert dimensionless.compute_reynolds_number(V, D, Nu_f) == pytest.approx(D) + + +@hypothesis.given(kinematic_viscosity=st.floats(min_value=1e-10, allow_nan=False)) +def test_reynolds_number_scales_inversly_with_kinematic_viscosity(kinematic_viscosity): + V = 1 + D = 1 + Nu_f = kinematic_viscosity + + assert dimensionless.compute_reynolds_number(V, D, Nu_f) == approx(1 / Nu_f) + + +def test_reynolds_number_with_example(): + V = 0.1 + D = 1.2 + Nu_f = 10 + + assert dimensionless.compute_reynolds_number(V, D, Nu_f) == approx(0.012) + + +@hypothesis.given(conductor_diameter=st.floats(min_value=1e-5, max_value=1e5, allow_nan=False)) +def test_grashof_number_scales_cubic_with_conductor_diameter(conductor_diameter): + D = conductor_diameter + T_s = 1 + T_a = 0 + Nu_f = 1 + g = 273.65 # 0.5*T_s in Kelvin + + Gr = dimensionless.compute_grashof_number(D, T_s, T_a, Nu_f, g) + assert Gr == approx(D**3) + + +@hypothesis.given(surface_temperature=st.floats(min_value=273, max_value=1000, allow_nan=False)) +def test_grashof_number_scales_correctly_with_surface_temperature(surface_temperature): + D = 1 + T_s = surface_temperature + T_a = 0 + Nu_f = 1 + g = 0.5 * T_s + 273.15 + + Gr = dimensionless.compute_grashof_number(D, T_s, T_a, Nu_f, g) + assert Gr == approx(abs(T_s)) + + +@hypothesis.given(air_temperature=st.floats(min_value=273, max_value=1000, allow_nan=False)) +def test_grashof_number_scales_correctly_with_air_temperature(air_temperature): + D = 1 + T_s = 0 + T_a = air_temperature + Nu_f = 1 + g = 0.5 * T_a + 273.15 + + Gr = dimensionless.compute_grashof_number(D, T_s, T_a, Nu_f, g) + assert Gr == approx(abs(T_a)) + + +@hypothesis.given( + kinematic_viscosity_of_air=st.floats(min_value=1e-5, max_value=1e10, allow_nan=False) +) +def test_grashof_number_scales_inversely_squared_with_kinematic_viscosity( + kinematic_viscosity_of_air, +): + D = 1 + T_s = 1 + T_a = 0 + Nu_f = kinematic_viscosity_of_air + g = 273.65 # 0.5*T_s in Kelvin + + Gr = dimensionless.compute_grashof_number(D, T_s, T_a, Nu_f, g) + assert Gr == approx(1 / (Nu_f**2)) + + +@hypothesis.given(coefficient_of_gravity=st.floats(min_value=1e-3, max_value=1e3, allow_nan=False)) +def test_grashof_number_scales_linearly_with_gravity(coefficient_of_gravity): + D = 1 + T_s = 1 + T_a = 0 + Nu_f = np.sqrt(1 / 273.65) + g = coefficient_of_gravity + + Gr = dimensionless.compute_grashof_number(D, T_s, T_a, Nu_f, g) + assert Gr == approx(g) + + +def test_grashof_number_with_example(): + D = 4 + T_s = 4 + T_a = 2 + Nu_f = 2 + g = 276.15 + + Gr = dimensionless.compute_grashof_number(D, T_s, T_a, Nu_f, g) + assert Gr == approx(32) + + +@hypothesis.given(thermal_conductivity=st.floats(min_value=1e-10, max_value=1e10, allow_nan=False)) +def test_prandtl_number_scales_inversely_with_thermal_conductivity(thermal_conductivity): + c_f = 1 + mu_f = 1 + lambda_f = thermal_conductivity + + assert dimensionless.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(1 / lambda_f) + + +@hypothesis.given(dynamic_viscosity=st.floats(min_value=1e-10, max_value=1e10, allow_nan=False)) +def test_prandtl_number_scales_linearly_with_dynamic_viscosity(dynamic_viscosity): + c_f = 1 + mu_f = dynamic_viscosity + lambda_f = 1 + + assert dimensionless.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(mu_f) + + +@hypothesis.given( + specific_heat_capacity=st.floats(min_value=1e-10, max_value=1e10, allow_nan=False) +) +def test_prandtl_number_scales_linearly_with_specific_heat_capacity(specific_heat_capacity): + c_f = specific_heat_capacity + mu_f = 1 + lambda_f = 1 + + assert dimensionless.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(c_f) + + +def test_compute_prandtl_number_with_example(): + c_f = 0.5 + mu_f = 3 + lambda_f = 0.5 + + assert dimensionless.compute_prandtl_number(lambda_f, mu_f, c_f) == approx(3) From 0580bc3c6484ec95d379d210dfee3f3aeb228527 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Thu, 30 May 2024 09:33:13 +0200 Subject: [PATCH 03/18] Complete implementation of Cigre 207, extract common solar heating --- .../equations/cigre207/convective_cooling.py | 1 - linerate/equations/cigre207/solar_heating.py | 99 ++---------- linerate/equations/cigre601/solar_heating.py | 34 +--- linerate/equations/solar_heating.py | 33 ++++ linerate/model.py | 4 +- linerate/models/Cigre207.md | 16 +- linerate/models/cigre207.py | 5 +- tests/equations/__init__.py | 0 .../cigre207/test_convective_cooling.py | 150 ++++++++++++++++++ .../equations/cigre601/test_solar_heating.py | 37 ----- tests/equations/test_convective_cooling.py | 36 +++++ tests/equations/test_radiative_cooling.py | 17 ++ tests/equations/test_solar_heating.py | 51 ++++++ 13 files changed, 318 insertions(+), 165 deletions(-) create mode 100644 linerate/equations/solar_heating.py create mode 100644 tests/equations/__init__.py create mode 100644 tests/equations/cigre207/test_convective_cooling.py create mode 100644 tests/equations/test_solar_heating.py diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py index 78e7215..8a0be0a 100644 --- a/linerate/equations/cigre207/convective_cooling.py +++ b/linerate/equations/cigre207/convective_cooling.py @@ -299,7 +299,6 @@ def _compute_horizontal_natural_nusselt_number( 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 diff --git a/linerate/equations/cigre207/solar_heating.py b/linerate/equations/cigre207/solar_heating.py index 5f0cf88..df756fa 100644 --- a/linerate/equations/cigre207/solar_heating.py +++ b/linerate/equations/cigre207/solar_heating.py @@ -1,39 +1,31 @@ import numpy as np from numpy import pi -from ...units import Meter, Unitless, WattPerMeter, WattPerSquareMeter +from .. import cigre601 +from ...units import Meter, Unitless, WattPerSquareMeter def compute_direct_solar_radiation( sin_solar_altitude: Unitless, - clearness_ratio: Unitless, height_above_sea_level: Meter, ) -> WattPerSquareMeter: r"""Compute the direct solar radiation. - Equation (10-11) on page 19 of :cite:p:`cigre601`. Equation (10) states that the direct solar + On page 38 of :cite:p:`cigre601`. Equation (10) states that the direct solar radiation on a surface normal to the solar beam at sea level, :math:`I_{B(0)}`, is given by .. math:: N_s \frac{1280 \sin(H_s)}{\sin(H_s) + 0.314}, - where :math:`N_s` is the clearness ratio which is used to adjust the amount of radiation - compared to what goes through a standard Indian atmosphere, and :math:`H_s` is the solar - altitude. - - While the solar radiation model is based on :cite:p:`sharma1965interrelationships` and - therefore have parameters estimated for an Indian atmosphere, it gives comparable results to - the solar radiation model in the IEEE standard :cite:p:`ieee738`. It is therefore reasonable to - assume that the parameters work in other climates as well. + where :math:`H_s` is the solar altitude. + To correct for height above sea level, we use the Eq. 19 from Cigre 601, since no equation is provided + in Cigre 207. Parameters ---------- sin_solar_altitude: :math:`\sin\left(H_s\right)`. The sine of the solar altitude. - clearness_ratio: - :math:`N_s`. The clearness ratio (or clearness number in - :cite:p:`sharma1965interrelationships,cigre207`). height_above_sea_level: :math:`y~\left[\text{m}\right]`. The conductor's altitude. @@ -41,40 +33,9 @@ def compute_direct_solar_radiation( ------- Union[float, float64, ndarray[Any, dtype[float64]]] :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. - - Note - ---- - The 1280 originates and 0.314 in the above equation originates from - :cite:p:`sharma1965interrelationships`, which is cited in :cite:p:`morgan1982thermal` (which is - listed as the reference in :cite:p:`cigre601`). In :cite:p:`sharma1965interrelationships` the - empirical relationship - - .. math:: - - I_{B(0)} = \frac{1.842 \sin(H_s)}{\sin(H_s) + 0.3135}~\text{Ly}~\text{min}^{-1} - - is introduced, and by converting from Langley per minute to :math:`\text{W}~\text{m}^{-2}`, we - obtain - - .. math:: - - I_{B(0)} = N_s \frac{1284.488 \sin(H_s)}{\sin(H_s) + 0.3135}~\text{W}~\text{m}^{-2}, - - which is equal to the equation we use (with three significant digits). """ - sin_H_s = sin_solar_altitude - N_s = clearness_ratio - y = height_above_sea_level - - I_B_0 = N_s * 1280 * sin_H_s / (sin_H_s + 0.314) - # Equation 19 says that - # I_B = I_B_0 * (1 + 1.4e-4 * y * (1367/I_B_0 - 1)) - # However, if I_B_0 = 0, this will divide by 0. To return NaN-values if and only - # if the input is NaN, we therefore reformulate it - scaled_y = 1.4e-4 * y - I_B = I_B_0 * (1 - scaled_y) + 1367 * scaled_y - - return np.where(sin_H_s >= 0, I_B, 0 * I_B) + clearness_ratio = 1.0 + return cigre601.solar_heating.compute_direct_solar_radiation(sin_solar_altitude, clearness_ratio, height_above_sea_level) def compute_diffuse_sky_radiation( @@ -83,9 +44,9 @@ def compute_diffuse_sky_radiation( ) -> WattPerSquareMeter: r"""Compute the diffuse radiation (light scattered in the atmosphere). - Equation (13) on page 20 of :cite:p:`cigre601`. + On page 38 of :cite:p:`cigre207`. - This equation differ from :cite:p:`cigre207`, however the difference is small, and the + This equation differ from :cite:p:`cigre601`, however the difference is small, and the diffuse radiation is a small contributor to the overall solar radiation, so the total discrepancy between the models is small. @@ -103,7 +64,7 @@ def compute_diffuse_sky_radiation( """ sin_H_s = sin_solar_altitude I_B = direct_solar_radiation - return np.maximum(0, (430.5 - 0.3288 * I_B)) * np.maximum(0, sin_H_s) + return np.maximum(0, (570 - 0.47 * I_B)) * np.maximum(0, sin_H_s)**1.2 def compute_global_radiation_intensity( @@ -115,14 +76,14 @@ def compute_global_radiation_intensity( ) -> WattPerSquareMeter: r"""Compute the global radiation intensity experienced by the conductor. - Equation (9) on page 18 of :cite:p:`cigre601` state that the global radiation intensity, + Equation (47) on page 38 of :cite:p:`cigre207` state that the global radiation intensity, :math:`I_T`, is given by .. math:: I_T = I_B \left(\sin(\eta) + 0.5 F \pi \sin(H_s)\right) + - I_d \left(1 + 0.5 F \pi\right), + 0.5 I_d \pi \left(1 + F \right), where :math:`\eta` is the incidence angle of the sun on the line, :math:`H_s` is the solar altitude and :math:`F` is the ground albedo (amount of radiation diffusely reflected from the @@ -181,36 +142,4 @@ def compute_global_radiation_intensity( sin_eta = sin_angle_of_sun_on_line F_pi_half = 0.5 * pi * F - return I_B * (sin_eta + F_pi_half * sin_H_s) + I_d * (1 + F_pi_half) # type: ignore - - -def compute_solar_heating( - absorptivity: Unitless, - global_radiation_intensity: WattPerSquareMeter, - conductor_diameter: Meter, -) -> WattPerMeter: - r"""Compute the solar heating experienced by the conductor. - - Equation (8) on page 18 of :cite:p:`cigre601`. - - Parameters - ---------- - absorptivity: - :math:`\alpha_s`. Material constant. According to :cite:p:`cigre601`, it starts at - approximately 0.2 for new cables and reaches a constant value of approximately 0.9 - after about one year. - global_radiation_intensity: - :math:`I_T~\left[\text{W}~\text{m}^{-2}\right]`.The global radiation intensity. - conductor_diameter: - :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_S~\left[\text{W}~\text{m}^{-1}\right]`. The solar heating of the conductor - """ - alpha_s = absorptivity - I_T = global_radiation_intensity - D = conductor_diameter - - return alpha_s * I_T * D + return I_B * (sin_eta + F_pi_half * sin_H_s) + I_d * pi/2 * (1 + F) # type: ignore diff --git a/linerate/equations/cigre601/solar_heating.py b/linerate/equations/cigre601/solar_heating.py index cb2c0bc..575920f 100644 --- a/linerate/equations/cigre601/solar_heating.py +++ b/linerate/equations/cigre601/solar_heating.py @@ -85,7 +85,7 @@ def compute_diffuse_sky_radiation( Equation (13) on page 20 of :cite:p:`cigre601`. - This equation differ from :cite:p:`cigre207`, however the difference is small, and the + This equation differs from :cite:p:`cigre207`, however the difference is small, and the diffuse radiation is a small contributor to the overall solar radiation, so the total discrepancy between the models is small. @@ -182,35 +182,3 @@ def compute_global_radiation_intensity( F_pi_half = 0.5 * pi * F return I_B * (sin_eta + F_pi_half * sin_H_s) + I_d * (1 + F_pi_half) # type: ignore - - -def compute_solar_heating( - absorptivity: Unitless, - global_radiation_intensity: WattPerSquareMeter, - conductor_diameter: Meter, -) -> WattPerMeter: - r"""Compute the solar heating experienced by the conductor. - - Equation (8) on page 18 of :cite:p:`cigre601`. - - Parameters - ---------- - absorptivity: - :math:`\alpha_s`. Material constant. According to :cite:p:`cigre601`, it starts at - approximately 0.2 for new cables and reaches a constant value of approximately 0.9 - after about one year. - global_radiation_intensity: - :math:`I_T~\left[\text{W}~\text{m}^{-2}\right]`.The global radiation intensity. - conductor_diameter: - :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_S~\left[\text{W}~\text{m}^{-1}\right]`. The solar heating of the conductor - """ - alpha_s = absorptivity - I_T = global_radiation_intensity - D = conductor_diameter - - return alpha_s * I_T * D diff --git a/linerate/equations/solar_heating.py b/linerate/equations/solar_heating.py new file mode 100644 index 0000000..0cfe465 --- /dev/null +++ b/linerate/equations/solar_heating.py @@ -0,0 +1,33 @@ +from linerate.units import WattPerSquareMeter, Unitless, Meter, WattPerMeter + + +def compute_solar_heating( + absorptivity: Unitless, + global_radiation_intensity: WattPerSquareMeter, + conductor_diameter: Meter, +) -> WattPerMeter: + r"""Compute the solar heating experienced by the conductor. + + Equation (8) on page 18 of :cite:p:`cigre601` and (11) on page 4 in :cite:p:`cigre207`. + + Parameters + ---------- + absorptivity: + :math:`\alpha_s`. Material constant. According to :cite:p:`cigre601`, it starts at + approximately 0.2 for new cables and reaches a constant value of approximately 0.9 + after about one year. + global_radiation_intensity: + :math:`I_T~\left[\text{W}~\text{m}^{-2}\right]`.The global radiation intensity. + conductor_diameter: + :math:`D~\left[\text{m}\right]`. Outer diameter of the conductor. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_S~\left[\text{W}~\text{m}^{-1}\right]`. The solar heating of the conductor + """ + alpha_s = absorptivity + I_T = global_radiation_intensity + D = conductor_diameter + + return alpha_s * I_T * D diff --git a/linerate/model.py b/linerate/model.py index 3543298..7d3af10 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -19,7 +19,7 @@ joule_heating, math, radiative_cooling, - solar_angles, convective_cooling, dimensionless, + solar_angles, convective_cooling, dimensionless, solar_heating, ) from linerate.equations.math import switch_cos_sin from linerate.types import Span, Weather @@ -355,7 +355,7 @@ def compute_solar_heating( I_T = cigre601.solar_heating.compute_global_radiation_intensity( I_B, I_d, F, sin_eta, sin_H_s ) - return cigre601.solar_heating.compute_solar_heating( + return solar_heating.compute_solar_heating( alpha_s, I_T, D, diff --git a/linerate/models/Cigre207.md b/linerate/models/Cigre207.md index 3912193..62c190a 100644 --- a/linerate/models/Cigre207.md +++ b/linerate/models/Cigre207.md @@ -15,16 +15,22 @@ We do not implement skin effect, since this is included in the AC resistance val ## Magnetic effect The effect of magnetic heating seems to be intertwined with the skin effect in Cigre 207. -To make things simpler and hopefully accurate enough, we use the same linearised model as our implementation of -Cigre 601. -**Maybe we have to come back to this one...** +To make things simpler and hopefully accurate enough, we use the same linearised model for magnetic effects in ACSR +lines as our implementation of Cigre 601. ## Solar heating -Solar heating is the same as for Cigre 601, but with clearness ratio 1. +Cigre 207 assumes that the diffuse radiation is uniformly directed, this differs from Cigre 601. +Apart from this, solar heating is the same as for Cigre 601, with clearness ratio 1. ## Convective cooling + There are a few inconsistencies in the Cigre 207 modelling of convective cooling. The definition of the Reynolds number is given as $Re = \rho_r V D/\nu_f$, where \rho_r is the air density relative to density at sea level. -This factor is not present in the ordinary definition of Reynolds number, so we choose to omit it. \ No newline at end of file +This factor is not present in the ordinary definition of Reynolds number, but seems to indicate that the kinematic +viscosity $\nu_f$ has to be corrected for the air relative density. + +Cigre 207 has a non-smooth transition to low wind speeds (<0.5m/s) at angles of attack less than 45 degrees. +This leads to an increase in cooling as the wind speed drops below 0.5m/s, which seems non-physical. +This discontinuity is removed in Cigre 601. \ No newline at end of file diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index 9f2ae92..5343364 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -1,7 +1,8 @@ from abc import abstractmethod from linerate import ThermalModel, Span, Weather -from linerate.equations import solar_angles, cigre601, math, cigre207, dimensionless, convective_cooling, joule_heating +from linerate.equations import solar_angles, cigre601, math, cigre207, dimensionless, convective_cooling, joule_heating, \ + solar_heating from linerate.equations.math import switch_cos_sin from linerate.model import _copy_method_docstring from linerate.units import Date, Celsius, Ampere, WattPerMeter, OhmPerMeter @@ -53,7 +54,7 @@ def compute_solar_heating( I_T = cigre601.solar_heating.compute_global_radiation_intensity( I_B, I_d, F, sin_eta, sin_H_s ) - return cigre601.solar_heating.compute_solar_heating( + return solar_heating.compute_solar_heating( alpha_s, I_T, D, diff --git a/tests/equations/__init__.py b/tests/equations/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/equations/cigre207/test_convective_cooling.py b/tests/equations/cigre207/test_convective_cooling.py new file mode 100644 index 0000000..e6b0a7f --- /dev/null +++ b/tests/equations/cigre207/test_convective_cooling.py @@ -0,0 +1,150 @@ +import numpy as np +from pytest import approx +from linerate.equations import cigre207, dimensionless, convective_cooling + + +def test_matches_example1(): + # See Appendix 1, Example 1 in Cigre 207 + y = 1600 + rho_r = cigre207.convective_cooling.compute_relative_air_density(y) + assert rho_r == approx(0.8306, rel=1e-4) + T_s = 57 + T_amb = 40 + T_f = (T_s + T_amb)/2 + nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) + assert nu_f == approx(1.78e-5, rel=1e-3) + v = 2 + D = 0.0286 + Re = cigre207.convective_cooling.compute_reynolds_number(v, D, nu_f, rho_r) + assert Re == approx(2670, rel=1e-3) + lambda_f = cigre207.convective_cooling.compute_thermal_conductivity_of_air(T_f) + assert lambda_f == approx(0.0277, 1e-3) + + +def test_matches_example1_nusselt_number(): + # See Appendix 1, Example 1 in Cigre 207 + D = 0.0286 + d = 0.00318 + Rs = dimensionless.compute_conductor_roughness(D, d) + Re = 2670 + Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) + assert Nu_90 == approx(26.45, 1e-4) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + assert Nu_45 == approx(22.34, rel=1e-4) + + +def test_matches_example2(): + # See Appendix 1, Example 2 in Cigre 207 + y = 1600 + rho_r = cigre207.convective_cooling.compute_relative_air_density(y) + T_s = 93 + T_amb = 40 + v = 0.2 + D = 0.0286 + T_f = (T_s + T_amb)/2 + nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) + nu_f_exp = 1.95e-5 + assert nu_f == approx(nu_f_exp, rel=1e-3) + Re = cigre207.convective_cooling.compute_reynolds_number(v, D, nu_f, rho_r) + assert Re == approx(243.8, rel=2e-3) + lambda_f = cigre207.convective_cooling.compute_thermal_conductivity_of_air(T_f) + assert lambda_f == approx(0.0290, 1e-3) + d = 0.00318 + Rs = dimensionless.compute_conductor_roughness(D, d) + Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) + assert Nu_90 == approx(8.53, 1e-3) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + assert Nu_45 == approx(7.20, rel=2e-3) + Nu_cor = cigre207.convective_cooling.compute_low_wind_speed_nusseltnumber(Nu_90) + assert Nu_cor == approx(4.69, rel=1e-3) + Gr = dimensionless.compute_grashof_number(D, T_s, T_amb, nu_f_exp) + assert Gr == approx(94387, rel=3e-3) + Pr = cigre207.convective_cooling.compute_prandtl_number(T_f) + assert Pr == approx(0.698, rel=1e-3) + Nu_natural = cigre207.convective_cooling.compute_horizontal_natural_nusselt_number(Gr, Pr) + assert Nu_natural == approx(7.69, rel=1e-3) + Nu_eff = cigre207.convective_cooling.compute_nusselt_number(Nu_45, Nu_natural, Nu_cor, v) + assert Nu_eff == approx(7.69, rel=1e-3) + + +def test_matches_example3(): + # See Appendix 1, Example 3 in Cigre 207 + y = 1600 + rho_r = cigre207.convective_cooling.compute_relative_air_density(y) + T_s = 75 + T_amb = 40 + v = 2 + D = 0.0286 + d = 0.00318 + T_f = (T_s + T_amb)/2 + nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) + nu_f_exp = 1.866e-5 + assert nu_f == approx(nu_f_exp, rel=1e-3) + Re = cigre207.convective_cooling.compute_reynolds_number(v, D, nu_f, rho_r) + assert Re == approx(2547.5, rel=1e-3) + lambda_f = cigre207.convective_cooling.compute_thermal_conductivity_of_air(T_f) + assert lambda_f == approx(0.0283, 2e-3) + Rs = dimensionless.compute_conductor_roughness(D, d) + Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) + assert Nu_90 == approx(25.77, 1e-3) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + assert Nu_45 == approx(21.8, rel=2e-3) + + +def test_matches_example4(): + # See Appendix 1, Example 4 in Cigre 207 + y = 1600 + rho_r = cigre207.convective_cooling.compute_relative_air_density(y) + T_s = 75 + T_amb = 40 + v = 0.4 + D = 0.0286 + T_f = (T_s + T_amb)/2 + nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) + nu_f_exp = 1.866e-5 + assert nu_f == approx(nu_f_exp, rel=1e-3) + Re = cigre207.convective_cooling.compute_reynolds_number(v, D, nu_f, rho_r) + assert Re == approx(509.6, rel=2e-3) + lambda_f = cigre207.convective_cooling.compute_thermal_conductivity_of_air(T_f) + assert lambda_f == approx(0.0283, 2e-3) + d = 0.00318 + Rs = dimensionless.compute_conductor_roughness(D, d) + Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) + assert Nu_90 == approx(12.08, 1e-3) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + assert Nu_45 == approx(10.2, rel=2e-3) + Nu_cor = cigre207.convective_cooling.compute_low_wind_speed_nusseltnumber(Nu_90) + assert Nu_cor == approx(6.64, rel=1e-3) + Gr = dimensionless.compute_grashof_number(D, T_s, T_amb, nu_f_exp) + assert Gr == approx(69922.7, rel=3e-3) + Pr = cigre207.convective_cooling.compute_prandtl_number(T_f) + assert Pr == approx(0.701, rel=1e-3) + Nu_natural = cigre207.convective_cooling.compute_horizontal_natural_nusselt_number(Gr, Pr) + assert Nu_natural == approx(7.14, rel=1e-3) + Nu_eff = cigre207.convective_cooling.compute_nusselt_number(Nu_45, Nu_natural, Nu_cor, v) + assert Nu_eff == approx(10.2, rel=1e-3) + + +def test_matches_example5(): + # See Appendix 1, Example 5 in Cigre 207 + y = 300 + rho_r = cigre207.convective_cooling.compute_relative_air_density(y) + assert rho_r == approx(0.966, rel=1e-3) + T_s = 57 + T_amb = 40 + v = 2 + D = 0.0286 + T_f = (T_s + T_amb)/2 + nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) + nu_f_exp = 1.78e-5 + assert nu_f == approx(nu_f_exp, rel=1e-3) + Re = cigre207.convective_cooling.compute_reynolds_number(v, D, nu_f, rho_r) + assert Re == approx(3106, rel=2e-3) + lambda_f = cigre207.convective_cooling.compute_thermal_conductivity_of_air(T_f) + assert lambda_f == approx(0.0277, 2e-3) + d = 0.00318 + Rs = dimensionless.compute_conductor_roughness(D, d) + Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) + assert Nu_90 == approx(29.85, 1e-3) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + assert Nu_45 == approx(25.21, rel=2e-3) diff --git a/tests/equations/cigre601/test_solar_heating.py b/tests/equations/cigre601/test_solar_heating.py index bd875fc..abe7091 100644 --- a/tests/equations/cigre601/test_solar_heating.py +++ b/tests/equations/cigre601/test_solar_heating.py @@ -218,40 +218,3 @@ def test_global_radiation_intensity_with_examples(): I_T = solar_heating.compute_global_radiation_intensity(I_B, I_d, F, sin_eta, sin_H_s) assert I_T == approx(2.5) - - -@hypothesis.given(conductor_diameter=st.floats(allow_nan=False)) -def test_solar_heating_scales_linearly_with_conductor_diameter(conductor_diameter): - D = conductor_diameter - alpha_s = 1 - I_T = 1 - - assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx(D) - - -@hypothesis.given(solar_absorptivity=st.floats(allow_nan=False)) -def test_solar_heating_scales_linearly_with_solar_absorptivity(solar_absorptivity): - D = 1 - alpha_s = solar_absorptivity - I_T = 1 - - assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx(alpha_s) - - -@hypothesis.given(global_radiation_intensity=st.floats(allow_nan=False)) -def test_solar_heating_scales_linearly_with_global_radiation_intensity(global_radiation_intensity): - D = 1 - alpha_s = 1 - I_T = global_radiation_intensity - - assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx( - global_radiation_intensity - ) - - -def test_solar_heating_scales_linearly_with_example(): - D = 0.6 - alpha_s = 0.5 - I_T = 2 - - assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx(0.6) diff --git a/tests/equations/test_convective_cooling.py b/tests/equations/test_convective_cooling.py index 876f749..f16f208 100644 --- a/tests/equations/test_convective_cooling.py +++ b/tests/equations/test_convective_cooling.py @@ -95,3 +95,39 @@ def test_convective_cooling_scales_affinely_with_surface_temperature(surface_tem lambda_f = 1 / np.pi assert convective_cooling.compute_convective_cooling(T_s, T_a, Nu, lambda_f) + 1 == approx(T_s) + + +def test_cooling_matches_cigre207_example1(): + # See Appendix 1, Example 1 in Cigre 207 + T_s = 57 + T_amb = 40 + Nu_45 = 22.34 + cooling = convective_cooling.compute_convective_cooling(T_s, T_amb, Nu_45, 0.0277) + assert cooling == approx(33.04, rel=1e-3) + + +def test_cooling_matches_cigre207_example2(): + # See Appendix 1, Example 2 in Cigre 207 + T_s = 93 + T_amb = 40 + Nu_45 = 7.69 + cooling = convective_cooling.compute_convective_cooling(T_s, T_amb, Nu_45, 0.0290) + assert cooling == approx(37.13, rel=1e-3) + + +def test_cooling_matches_cigre207_example3(): + # See Appendix 1, Example 3 in Cigre 207 + T_s = 75 + T_amb = 40 + Nu_45 = 21.8 + cooling = convective_cooling.compute_convective_cooling(T_s, T_amb, Nu_45, 0.0283) + assert cooling == approx(67.8, rel=1e-3) + + +def test_cooling_matches_cigre207_example4(): + # See Appendix 1, Example 4 in Cigre 207 + T_s = 75 + T_amb = 40 + Nu_45 = 10.2 + cooling = convective_cooling.compute_convective_cooling(T_s, T_amb, Nu_45, 0.0283) + assert cooling == approx(31.78, rel=2e-3) \ No newline at end of file diff --git a/tests/equations/test_radiative_cooling.py b/tests/equations/test_radiative_cooling.py index 9b8211d..6170138 100644 --- a/tests/equations/test_radiative_cooling.py +++ b/tests/equations/test_radiative_cooling.py @@ -83,3 +83,20 @@ def test_radiative_cooling_with_example( T_a = air_temperature assert radiative_cooling.compute_radiative_cooling(T_s, T_a, D, epsilon) == approx(cooling) + + +@pytest.mark.parametrize( + "surface_temperature, expected_cooling", + [ + (57, 5.76), + (93, 21.27), + (75, 12.92) + ], +) +def test_cooling_matches_cigre207_examples(surface_temperature, expected_cooling): + # See Appendix 1, Example 1 in Cigre 207 + T_a = 40 + D = 0.0286 + epsilon = 0.5 + cooling = radiative_cooling.compute_radiative_cooling(surface_temperature, T_a, D, epsilon) + assert cooling == approx(expected_cooling, rel=2e-3) \ No newline at end of file diff --git a/tests/equations/test_solar_heating.py b/tests/equations/test_solar_heating.py new file mode 100644 index 0000000..864366d --- /dev/null +++ b/tests/equations/test_solar_heating.py @@ -0,0 +1,51 @@ +import hypothesis +import hypothesis.strategies as st +from pytest import approx + +from linerate.equations import solar_heating + + +def test_matches_solar_heating_example_from_cigre207(): + # See Appendix 1, Example 1 in Cigre 207 + D = 0.0286 + global_solar_radiation = 980 + absorptivity = 0.5 + heating = solar_heating.compute_solar_heating(absorptivity, global_solar_radiation, D) + assert heating == approx(14.02, rel=1e-3) + + +@hypothesis.given(conductor_diameter=st.floats(allow_nan=False)) +def test_solar_heating_scales_linearly_with_conductor_diameter(conductor_diameter): + D = conductor_diameter + alpha_s = 1 + I_T = 1 + + assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx(D) + + +@hypothesis.given(solar_absorptivity=st.floats(allow_nan=False)) +def test_solar_heating_scales_linearly_with_solar_absorptivity(solar_absorptivity): + D = 1 + alpha_s = solar_absorptivity + I_T = 1 + + assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx(alpha_s) + + +@hypothesis.given(global_radiation_intensity=st.floats(allow_nan=False)) +def test_solar_heating_scales_linearly_with_global_radiation_intensity(global_radiation_intensity): + D = 1 + alpha_s = 1 + I_T = global_radiation_intensity + + assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx( + global_radiation_intensity + ) + + +def test_solar_heating_scales_linearly_with_example(): + D = 0.6 + alpha_s = 0.5 + I_T = 2 + + assert solar_heating.compute_solar_heating(alpha_s, I_T, D) == approx(0.6) From 0906eff02782f3c3ec30eccc849fd11234096f92 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Thu, 30 May 2024 10:06:37 +0200 Subject: [PATCH 04/18] Fix formatting, import Cigre207 in model.py --- linerate/equations/cigre207/__init__.py | 2 +- linerate/equations/cigre207/ac_resistance.py | 4 +- .../equations/cigre207/convective_cooling.py | 29 +++--- linerate/equations/cigre207/solar_heating.py | 14 +-- .../equations/cigre601/convective_cooling.py | 93 +++++++++---------- linerate/equations/cigre601/solar_heating.py | 2 +- linerate/equations/convective_cooling.py | 15 +-- linerate/equations/dimensionless.py | 20 +++- linerate/equations/solar_heating.py | 2 +- linerate/model.py | 8 +- linerate/models/Cigre207.md | 2 +- linerate/models/cigre207.py | 25 +++-- .../cigre207/test_convective_cooling.py | 33 ++++--- .../cigre601/test_convective_cooling.py | 3 - tests/equations/test_convective_cooling.py | 4 +- tests/equations/test_radiative_cooling.py | 8 +- 16 files changed, 142 insertions(+), 122 deletions(-) diff --git a/linerate/equations/cigre207/__init__.py b/linerate/equations/cigre207/__init__.py index af02088..d2c20ed 100644 --- a/linerate/equations/cigre207/__init__.py +++ b/linerate/equations/cigre207/__init__.py @@ -2,4 +2,4 @@ This submodule contains implementations of equations listed in :cite:p:`cigre207`. """ -from . import convective_cooling, solar_heating # noqa \ No newline at end of file +from . import convective_cooling, solar_heating # noqa diff --git a/linerate/equations/cigre207/ac_resistance.py b/linerate/equations/cigre207/ac_resistance.py index 748baef..550128f 100644 --- a/linerate/equations/cigre207/ac_resistance.py +++ b/linerate/equations/cigre207/ac_resistance.py @@ -1,5 +1,3 @@ -from typing import Union - from linerate.units import OhmPerMeter @@ -24,4 +22,4 @@ def correct_resistance_for_skin_effect( return 1.0123 * dc_resistance -# TODO: Implement section 2.1.2? \ No newline at end of file +# TODO: Implement section 2.1.2? diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py index 8a0be0a..3a727d4 100644 --- a/linerate/equations/cigre207/convective_cooling.py +++ b/linerate/equations/cigre207/convective_cooling.py @@ -9,11 +9,11 @@ Meter, MeterPerSecond, Radian, + SquareMeterPerSecond, Unitless, - WattPerMeterPerKelvin, SquareMeterPerSecond, + WattPerMeterPerKelvin, ) - # Physical quantities ##################### @@ -40,10 +40,9 @@ def compute_thermal_conductivity_of_air(film_temperature: Celsius) -> WattPerMet 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. +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`. @@ -58,7 +57,7 @@ def compute_relative_air_density( :math:`\rho_r`. The relative mass density of air. """ y = height_above_sea_level - return np.exp(-1.16e-4*y) + return np.exp(-1.16e-4 * y) def compute_kinematic_viscosity_of_air(film_temperature: Celsius) -> KilogramPerCubeMeter: @@ -79,8 +78,7 @@ def compute_kinematic_viscosity_of_air(film_temperature: Celsius) -> KilogramPer :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 - + return 1.32e-5 + 9.5e-8 * T_f def compute_prandtl_number( @@ -104,14 +102,14 @@ def compute_prandtl_number( Union[float, float64, ndarray[Any, dtype[float64]]] :math:`\text{Pr}`. The Prandtl number. """ - return 0.715 - 2.5e-4*film_temperature + return 0.715 - 2.5e-4 * film_temperature def compute_reynolds_number( wind_speed: MeterPerSecond, conductor_diameter: Meter, kinematic_viscosity_of_air: SquareMeterPerSecond, - relative_air_density: Unitless + relative_air_density: Unitless, ) -> Unitless: r"""Compute the Reynolds number using the conductor diameter as characteristic length scale. @@ -142,7 +140,6 @@ def compute_reynolds_number( return rho_r * v * D / nu_f - ## Nusselt number calculation ############################# @@ -219,8 +216,7 @@ def compute_low_wind_speed_nusseltnumber( 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 - + return 0.55 * perpendicular_flow_nusselt_number @vectorize(nopython=True) @@ -265,8 +261,7 @@ def correct_wind_direction_effect_on_nusselt_number( :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 + perpendicular_flow_nusselt_number, angle_of_attack ) @@ -337,7 +332,7 @@ def compute_nusselt_number( forced_convection_nusselt_number: Unitless, natural_nusselt_number: Unitless, low_wind_nusselt_number: Unitless, - wind_speed: MeterPerSecond + wind_speed: MeterPerSecond, ) -> Unitless: r"""Compute the nusselt number. diff --git a/linerate/equations/cigre207/solar_heating.py b/linerate/equations/cigre207/solar_heating.py index df756fa..ef6b4d0 100644 --- a/linerate/equations/cigre207/solar_heating.py +++ b/linerate/equations/cigre207/solar_heating.py @@ -1,8 +1,8 @@ import numpy as np from numpy import pi -from .. import cigre601 from ...units import Meter, Unitless, WattPerSquareMeter +from .. import cigre601 def compute_direct_solar_radiation( @@ -19,8 +19,8 @@ def compute_direct_solar_radiation( N_s \frac{1280 \sin(H_s)}{\sin(H_s) + 0.314}, where :math:`H_s` is the solar altitude. - To correct for height above sea level, we use the Eq. 19 from Cigre 601, since no equation is provided - in Cigre 207. + To correct for height above sea level, we use the Eq. 19 from Cigre 601, + since no equation is provided in Cigre 207. Parameters ---------- @@ -35,7 +35,9 @@ def compute_direct_solar_radiation( :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. """ clearness_ratio = 1.0 - return cigre601.solar_heating.compute_direct_solar_radiation(sin_solar_altitude, clearness_ratio, height_above_sea_level) + return cigre601.solar_heating.compute_direct_solar_radiation( + sin_solar_altitude, clearness_ratio, height_above_sea_level + ) def compute_diffuse_sky_radiation( @@ -64,7 +66,7 @@ def compute_diffuse_sky_radiation( """ sin_H_s = sin_solar_altitude I_B = direct_solar_radiation - return np.maximum(0, (570 - 0.47 * I_B)) * np.maximum(0, sin_H_s)**1.2 + return np.maximum(0, (570 - 0.47 * I_B)) * np.maximum(0, sin_H_s) ** 1.2 def compute_global_radiation_intensity( @@ -142,4 +144,4 @@ def compute_global_radiation_intensity( sin_eta = sin_angle_of_sun_on_line F_pi_half = 0.5 * pi * F - return I_B * (sin_eta + F_pi_half * sin_H_s) + I_d * pi/2 * (1 + F) # type: ignore + return I_B * (sin_eta + F_pi_half * sin_H_s) + I_d * pi / 2 * (1 + F) # type: ignore diff --git a/linerate/equations/cigre601/convective_cooling.py b/linerate/equations/cigre601/convective_cooling.py index c02ff29..2069062 100644 --- a/linerate/equations/cigre601/convective_cooling.py +++ b/linerate/equations/cigre601/convective_cooling.py @@ -15,16 +15,15 @@ WattPerMeterPerKelvin, ) - # Physical quantities ##################### def compute_temperature_gradient( - total_heat_gain: WattPerMeter, - conductor_thermal_conductivity: WattPerMeterPerKelvin, - core_diameter: Meter, - conductor_diameter: Meter, + total_heat_gain: WattPerMeter, + conductor_thermal_conductivity: WattPerMeterPerKelvin, + core_diameter: Meter, + conductor_diameter: Meter, ) -> Celsius: r"""Compute the difference between the core and surface temperature. @@ -64,8 +63,8 @@ def compute_temperature_gradient( if D_1 == 0: # TODO: Maybe lower tolerance? return 0.5 * tmp else: - D_1_sq = D_1 ** 2 - delta_D_sq = D ** 2 - D_1_sq + D_1_sq = D_1**2 + delta_D_sq = D**2 - D_1_sq return tmp * (0.5 - (D_1_sq / delta_D_sq) * np.log(D / D_1)) @@ -90,11 +89,11 @@ def compute_thermal_conductivity_of_air(film_temperature: Celsius) -> WattPerMet conductivity of air at the given temperature. """ T_f = film_temperature - return 2.368e-2 + 7.23e-5 * T_f - 2.763e-8 * (T_f ** 2) + return 2.368e-2 + 7.23e-5 * T_f - 2.763e-8 * (T_f**2) def compute_air_density( - film_temperature: Celsius, height_above_sea_level: Meter + film_temperature: Celsius, height_above_sea_level: Meter ) -> KilogramPerCubeMeter: r"""Approximation of the density of air at a given temperature and altitude. @@ -116,7 +115,7 @@ def compute_air_density( """ T_f = film_temperature y = height_above_sea_level - return (1.293 - 1.525e-4 * y + 6.379e-9 * (y ** 2)) / (1 + 0.00367 * T_f) + return (1.293 - 1.525e-4 * y + 6.379e-9 * (y**2)) / (1 + 0.00367 * T_f) def compute_dynamic_viscosity_of_air(film_temperature: Celsius) -> KilogramPerMeterPerSecond: @@ -138,11 +137,11 @@ def compute_dynamic_viscosity_of_air(film_temperature: Celsius) -> KilogramPerMe of air. """ T_f = film_temperature - return 17.239e-6 + 4.635e-8 * T_f - 2.03e-11 * (T_f ** 2) + return 17.239e-6 + 4.635e-8 * T_f - 2.03e-11 * (T_f**2) def compute_kinematic_viscosity_of_air( - dynamic_viscosity_of_air: KilogramPerMeterPerSecond, air_density: KilogramPerCubeMeter + dynamic_viscosity_of_air: KilogramPerMeterPerSecond, air_density: KilogramPerCubeMeter ) -> SquareMeterPerSecond: r"""Compute the kinematic viscosity of air. @@ -185,8 +184,8 @@ def _check_perpendicular_flow_nusseltnumber_out_of_bounds(reynolds_number, condu @vectorize(nopython=True) def _compute_perpendicular_flow_nusseltnumber( - reynolds_number: Unitless, - conductor_roughness: Meter, + reynolds_number: Unitless, + conductor_roughness: Meter, ) -> Unitless: # TODO: Look at references for this table Re = reynolds_number @@ -216,12 +215,12 @@ def _compute_perpendicular_flow_nusseltnumber( else: B, n = 0.048, 0.800 - return B * Re ** n # type: ignore + 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. @@ -253,9 +252,9 @@ def compute_perpendicular_flow_nusseltnumber( @vectorize(nopython=True) def _correct_wind_direction_effect_on_nusselt_number( - perpendicular_flow_nusselt_number: Unitless, - angle_of_attack: Radian, - conductor_roughness: Unitless, + perpendicular_flow_nusselt_number: Unitless, + angle_of_attack: Radian, + conductor_roughness: Unitless, ) -> Unitless: delta = angle_of_attack Nu_90 = perpendicular_flow_nusselt_number @@ -264,23 +263,23 @@ def _correct_wind_direction_effect_on_nusselt_number( sin_delta = np.sin(delta) if Rs == 0 or np.isnan(Rs): - sin_delta_sq = sin_delta ** 2 + sin_delta_sq = sin_delta**2 cos_delta_sq = 1 - sin_delta_sq correction_factor = (sin_delta_sq + cos_delta_sq * 0.0169) ** 0.225 else: if delta <= np.radians(24): - correction_factor = 0.42 + 0.68 * (sin_delta ** 1.08) + correction_factor = 0.42 + 0.68 * (sin_delta**1.08) else: - correction_factor = 0.42 + 0.58 * (sin_delta ** 0.90) + 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, - conductor_roughness: Unitless, + perpendicular_flow_nusselt_number: Unitless, + angle_of_attack: Radian, + conductor_roughness: Unitless, ) -> Unitless: r"""Correct the Nusselt number for the wind's angle-of-attack. @@ -315,7 +314,7 @@ def correct_wind_direction_effect_on_nusselt_number( 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): @@ -326,26 +325,26 @@ def _check_horizontal_natural_nusselt_number( @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 < 1e-1: return 0 elif GrPr < 1e2: - return 1.020 * GrPr ** 0.148 + return 1.020 * GrPr**0.148 elif GrPr < 1e4: - return 0.850 * GrPr ** 0.188 + return 0.850 * GrPr**0.188 elif GrPr < 1e7: - return 0.480 * GrPr ** 0.250 + return 0.480 * GrPr**0.250 else: - return 0.125 * GrPr ** 0.333 + return 0.125 * GrPr**0.333 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. @@ -378,8 +377,8 @@ def compute_horizontal_natural_nusselt_number( def _check_conductor_inclination( - conductor_inclination: Radian, - conductor_roughness: Unitless, + conductor_inclination: Radian, + conductor_roughness: Unitless, ) -> None: beta = np.degrees(conductor_inclination) Rs = conductor_roughness @@ -396,24 +395,24 @@ def _check_conductor_inclination( @vectorize(nopython=True) def _correct_natural_nusselt_number_inclination( - horizontal_natural_nusselt_number: Unitless, - conductor_inclination: Radian, - conductor_roughness: Unitless, + horizontal_natural_nusselt_number: Unitless, + conductor_inclination: Radian, + conductor_roughness: Unitless, ) -> Unitless: beta = np.degrees(conductor_inclination) Nu_nat = horizontal_natural_nusselt_number Rs = conductor_roughness if Rs == 0 or np.isnan(Rs): - return Nu_nat * (1 - 1.58e-4 * beta ** 1.5) + return Nu_nat * (1 - 1.58e-4 * beta**1.5) else: - return Nu_nat * (1 - 1.76e-6 * beta ** 2.5) + return Nu_nat * (1 - 1.76e-6 * beta**2.5) def correct_natural_nusselt_number_inclination( - horizontal_natural_nusselt_number: Unitless, - conductor_inclination: Radian, - conductor_roughness: Unitless, + horizontal_natural_nusselt_number: Unitless, + conductor_inclination: Radian, + conductor_roughness: Unitless, ) -> Unitless: r"""Correct the natural Nusselt number for the effect of the span inclination. @@ -447,8 +446,8 @@ def correct_natural_nusselt_number_inclination( def compute_nusselt_number( - forced_convection_nusselt_number: Unitless, - natural_nusselt_number: Unitless, + forced_convection_nusselt_number: Unitless, + natural_nusselt_number: Unitless, ) -> Unitless: r"""Compute the nusselt number. diff --git a/linerate/equations/cigre601/solar_heating.py b/linerate/equations/cigre601/solar_heating.py index 575920f..aef33c8 100644 --- a/linerate/equations/cigre601/solar_heating.py +++ b/linerate/equations/cigre601/solar_heating.py @@ -1,7 +1,7 @@ import numpy as np from numpy import pi -from ...units import Meter, Unitless, WattPerMeter, WattPerSquareMeter +from ...units import Meter, Unitless, WattPerSquareMeter def compute_direct_solar_radiation( diff --git a/linerate/equations/convective_cooling.py b/linerate/equations/convective_cooling.py index cec9a54..2c7b427 100644 --- a/linerate/equations/convective_cooling.py +++ b/linerate/equations/convective_cooling.py @@ -1,17 +1,18 @@ import numpy as np -from linerate.units import Celsius, Unitless, WattPerMeterPerKelvin, WattPerMeter +from linerate.units import Celsius, Unitless, WattPerMeter, WattPerMeterPerKelvin def compute_convective_cooling( - surface_temperature: Celsius, - air_temperature: Celsius, - nusselt_number: Unitless, - thermal_conductivity_of_air: WattPerMeterPerKelvin, + surface_temperature: Celsius, + air_temperature: Celsius, + nusselt_number: Unitless, + thermal_conductivity_of_air: WattPerMeterPerKelvin, ) -> WattPerMeter: r"""Compute the convective cooling of the conductor. - Equation (17) on page 24 of :cite:p:`cigre601` and Equation (12) on page 6 of :cite:p:`cigre207`. + Equation (17) on page 24 of :cite:p:`cigre601` + and Equation (12) on page 6 of :cite:p:`cigre207`. Parameters ---------- @@ -37,4 +38,4 @@ def compute_convective_cooling( T_a = air_temperature Nu = nusselt_number - return pi * lambda_f * (T_s - T_a) * Nu \ No newline at end of file + return pi * lambda_f * (T_s - T_a) * Nu diff --git a/linerate/equations/dimensionless.py b/linerate/equations/dimensionless.py index 10430ad..a1e7c3e 100644 --- a/linerate/equations/dimensionless.py +++ b/linerate/equations/dimensionless.py @@ -2,8 +2,17 @@ import numpy as np -from linerate.units import Unitless, WattPerMeterPerKelvin, KilogramPerMeterPerSecond, JoulePerKilogramPerKelvin, \ - MeterPerSecond, Meter, SquareMeterPerSecond, Celsius, MeterPerSquareSecond +from linerate.units import ( + Celsius, + JoulePerKilogramPerKelvin, + KilogramPerMeterPerSecond, + Meter, + MeterPerSecond, + MeterPerSquareSecond, + SquareMeterPerSecond, + Unitless, + WattPerMeterPerKelvin, +) def compute_reynolds_number( @@ -13,7 +22,8 @@ def compute_reynolds_number( ) -> Unitless: r"""Compute the Reynolds number using the conductor diameter as characteristic length scale. - Defined in the text on page 25 of :cite:p:`cigre601` and equation (2c) on page 10 in :cite:p:`ieee738`. + Defined in the text on page 25 of :cite:p:`cigre601` and equation (2c) + on page 10 in :cite:p:`ieee738`. The Reynolds number is a dimensionless quantity that can be used to assess if a stream is likely to be turbulent or not. It is given by @@ -125,8 +135,8 @@ def compute_prandtl_number( def compute_conductor_roughness( - conductor_diameter: Meter, - outer_layer_strand_diameter: Meter, + conductor_diameter: Meter, + outer_layer_strand_diameter: Meter, ) -> Unitless: r"""Compute the surface roughness of the conductor. diff --git a/linerate/equations/solar_heating.py b/linerate/equations/solar_heating.py index 0cfe465..a177116 100644 --- a/linerate/equations/solar_heating.py +++ b/linerate/equations/solar_heating.py @@ -1,4 +1,4 @@ -from linerate.units import WattPerSquareMeter, Unitless, Meter, WattPerMeter +from linerate.units import Meter, Unitless, WattPerMeter, WattPerSquareMeter def compute_solar_heating( diff --git a/linerate/model.py b/linerate/model.py index 7d3af10..9699124 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -11,15 +11,19 @@ from typing import Dict import numpy as np +from models.cigre207 import Cigre207 from linerate import solver from linerate.equations import ( cigre601, + convective_cooling, + dimensionless, ieee738, joule_heating, math, radiative_cooling, - solar_angles, convective_cooling, dimensionless, solar_heating, + solar_angles, + solar_heating, ) from linerate.equations.math import switch_cos_sin from linerate.types import Span, Weather @@ -32,7 +36,7 @@ WattPerMeter, ) -__all__ = ["ThermalModel", "Cigre601", "IEEE738"] +__all__ = ["ThermalModel", "Cigre601", "IEEE738", "Cigre207"] def _copy_method_docstring(parent_class): diff --git a/linerate/models/Cigre207.md b/linerate/models/Cigre207.md index 62c190a..ad23e73 100644 --- a/linerate/models/Cigre207.md +++ b/linerate/models/Cigre207.md @@ -33,4 +33,4 @@ viscosity $\nu_f$ has to be corrected for the air relative density. Cigre 207 has a non-smooth transition to low wind speeds (<0.5m/s) at angles of attack less than 45 degrees. This leads to an increase in cooling as the wind speed drops below 0.5m/s, which seems non-physical. -This discontinuity is removed in Cigre 601. \ No newline at end of file +This discontinuity is removed in Cigre 601. diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index 5343364..96277f2 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -1,11 +1,19 @@ from abc import abstractmethod -from linerate import ThermalModel, Span, Weather -from linerate.equations import solar_angles, cigre601, math, cigre207, dimensionless, convective_cooling, joule_heating, \ - solar_heating +from linerate import Span, ThermalModel, Weather +from linerate.equations import ( + cigre207, + cigre601, + convective_cooling, + dimensionless, + joule_heating, + math, + solar_angles, + solar_heating, +) from linerate.equations.math import switch_cos_sin from linerate.model import _copy_method_docstring -from linerate.units import Date, Celsius, Ampere, WattPerMeter, OhmPerMeter +from linerate.units import Ampere, Celsius, Date, OhmPerMeter, WattPerMeter class Cigre207(ThermalModel): @@ -80,9 +88,6 @@ def compute_convective_cooling( ) # Compute unitless quantities - # Reynolds number is defined in the text on page 5 of :cite:p:`cigre207`. - # The definition includes a relative air density, which does not make sense, so we omit it here and use the - # standard definition of Reynolds number instead. rho_r = cigre207.convective_cooling.compute_relative_air_density(y) Re = cigre207.convective_cooling.compute_reynolds_number(V, D, nu_f, rho_r) Gr = dimensionless.compute_grashof_number(D, T_c, T_a, nu_f) @@ -101,8 +106,10 @@ def compute_convective_cooling( Nu_0 = cigre207.convective_cooling.compute_horizontal_natural_nusselt_number(Gr, Pr) Nu = cigre207.convective_cooling.compute_nusselt_number( - forced_convection_nusselt_number=Nu_delta, natural_nusselt_number=Nu_0, - low_wind_nusselt_number=Nu_cor, wind_speed=V + forced_convection_nusselt_number=Nu_delta, + natural_nusselt_number=Nu_0, + low_wind_nusselt_number=Nu_cor, + wind_speed=V, ) return convective_cooling.compute_convective_cooling( diff --git a/tests/equations/cigre207/test_convective_cooling.py b/tests/equations/cigre207/test_convective_cooling.py index e6b0a7f..44695fb 100644 --- a/tests/equations/cigre207/test_convective_cooling.py +++ b/tests/equations/cigre207/test_convective_cooling.py @@ -1,6 +1,7 @@ import numpy as np from pytest import approx -from linerate.equations import cigre207, dimensionless, convective_cooling + +from linerate.equations import cigre207, dimensionless def test_matches_example1(): @@ -10,7 +11,7 @@ def test_matches_example1(): assert rho_r == approx(0.8306, rel=1e-4) T_s = 57 T_amb = 40 - T_f = (T_s + T_amb)/2 + T_f = (T_s + T_amb) / 2 nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) assert nu_f == approx(1.78e-5, rel=1e-3) v = 2 @@ -29,7 +30,9 @@ def test_matches_example1_nusselt_number(): Re = 2670 Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) assert Nu_90 == approx(26.45, 1e-4) - Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number( + Nu_90, np.radians(45) + ) assert Nu_45 == approx(22.34, rel=1e-4) @@ -41,7 +44,7 @@ def test_matches_example2(): T_amb = 40 v = 0.2 D = 0.0286 - T_f = (T_s + T_amb)/2 + T_f = (T_s + T_amb) / 2 nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) nu_f_exp = 1.95e-5 assert nu_f == approx(nu_f_exp, rel=1e-3) @@ -53,7 +56,9 @@ def test_matches_example2(): Rs = dimensionless.compute_conductor_roughness(D, d) Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) assert Nu_90 == approx(8.53, 1e-3) - Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number( + Nu_90, np.radians(45) + ) assert Nu_45 == approx(7.20, rel=2e-3) Nu_cor = cigre207.convective_cooling.compute_low_wind_speed_nusseltnumber(Nu_90) assert Nu_cor == approx(4.69, rel=1e-3) @@ -76,7 +81,7 @@ def test_matches_example3(): v = 2 D = 0.0286 d = 0.00318 - T_f = (T_s + T_amb)/2 + T_f = (T_s + T_amb) / 2 nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) nu_f_exp = 1.866e-5 assert nu_f == approx(nu_f_exp, rel=1e-3) @@ -87,7 +92,9 @@ def test_matches_example3(): Rs = dimensionless.compute_conductor_roughness(D, d) Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) assert Nu_90 == approx(25.77, 1e-3) - Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number( + Nu_90, np.radians(45) + ) assert Nu_45 == approx(21.8, rel=2e-3) @@ -99,7 +106,7 @@ def test_matches_example4(): T_amb = 40 v = 0.4 D = 0.0286 - T_f = (T_s + T_amb)/2 + T_f = (T_s + T_amb) / 2 nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) nu_f_exp = 1.866e-5 assert nu_f == approx(nu_f_exp, rel=1e-3) @@ -111,7 +118,9 @@ def test_matches_example4(): Rs = dimensionless.compute_conductor_roughness(D, d) Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) assert Nu_90 == approx(12.08, 1e-3) - Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number( + Nu_90, np.radians(45) + ) assert Nu_45 == approx(10.2, rel=2e-3) Nu_cor = cigre207.convective_cooling.compute_low_wind_speed_nusseltnumber(Nu_90) assert Nu_cor == approx(6.64, rel=1e-3) @@ -134,7 +143,7 @@ def test_matches_example5(): T_amb = 40 v = 2 D = 0.0286 - T_f = (T_s + T_amb)/2 + T_f = (T_s + T_amb) / 2 nu_f = cigre207.convective_cooling.compute_kinematic_viscosity_of_air(T_f) nu_f_exp = 1.78e-5 assert nu_f == approx(nu_f_exp, rel=1e-3) @@ -146,5 +155,7 @@ def test_matches_example5(): Rs = dimensionless.compute_conductor_roughness(D, d) Nu_90 = cigre207.convective_cooling.compute_perpendicular_flow_nusseltnumber(Re, Rs) assert Nu_90 == approx(29.85, 1e-3) - Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number(Nu_90, np.radians(45)) + Nu_45 = cigre207.convective_cooling.correct_wind_direction_effect_on_nusselt_number( + Nu_90, np.radians(45) + ) assert Nu_45 == approx(25.21, rel=2e-3) diff --git a/tests/equations/cigre601/test_convective_cooling.py b/tests/equations/cigre601/test_convective_cooling.py index 6e690c4..841f0f8 100644 --- a/tests/equations/cigre601/test_convective_cooling.py +++ b/tests/equations/cigre601/test_convective_cooling.py @@ -7,7 +7,6 @@ import linerate.equations.cigre601.convective_cooling as convective_cooling - # Tests for physical quantities ############################### @@ -260,8 +259,6 @@ def test_kinematic_viscosity_of_air_with_example(): assert convective_cooling.compute_kinematic_viscosity_of_air(mu_f, gamma) == approx(1) - - ## Nusselt number calculation ############################# diff --git a/tests/equations/test_convective_cooling.py b/tests/equations/test_convective_cooling.py index f16f208..b0fc0dc 100644 --- a/tests/equations/test_convective_cooling.py +++ b/tests/equations/test_convective_cooling.py @@ -1,8 +1,8 @@ import hypothesis import hypothesis.strategies as st import numpy as np - from pytest import approx + from linerate.equations import convective_cooling @@ -130,4 +130,4 @@ def test_cooling_matches_cigre207_example4(): T_amb = 40 Nu_45 = 10.2 cooling = convective_cooling.compute_convective_cooling(T_s, T_amb, Nu_45, 0.0283) - assert cooling == approx(31.78, rel=2e-3) \ No newline at end of file + assert cooling == approx(31.78, rel=2e-3) diff --git a/tests/equations/test_radiative_cooling.py b/tests/equations/test_radiative_cooling.py index 6170138..27586bf 100644 --- a/tests/equations/test_radiative_cooling.py +++ b/tests/equations/test_radiative_cooling.py @@ -87,11 +87,7 @@ def test_radiative_cooling_with_example( @pytest.mark.parametrize( "surface_temperature, expected_cooling", - [ - (57, 5.76), - (93, 21.27), - (75, 12.92) - ], + [(57, 5.76), (93, 21.27), (75, 12.92)], ) def test_cooling_matches_cigre207_examples(surface_temperature, expected_cooling): # See Appendix 1, Example 1 in Cigre 207 @@ -99,4 +95,4 @@ def test_cooling_matches_cigre207_examples(surface_temperature, expected_cooling D = 0.0286 epsilon = 0.5 cooling = radiative_cooling.compute_radiative_cooling(surface_temperature, T_a, D, epsilon) - assert cooling == approx(expected_cooling, rel=2e-3) \ No newline at end of file + assert cooling == approx(expected_cooling, rel=2e-3) From 39d1bbd320801bdca0e5b54ae7e145a001fc88c5 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Thu, 30 May 2024 13:30:22 +0200 Subject: [PATCH 05/18] Fix circular imports --- linerate/model.py | 312 +++---------------------------- linerate/models/cigre207.py | 4 +- linerate/models/thermal_model.py | 281 ++++++++++++++++++++++++++++ 3 files changed, 304 insertions(+), 293 deletions(-) create mode 100644 linerate/models/thermal_model.py diff --git a/linerate/model.py b/linerate/model.py index 9699124..9c3aa35 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -6,26 +6,22 @@ ``linerate.solver`` modules. """ -from abc import ABC, abstractmethod from numbers import Real -from typing import Dict import numpy as np -from models.cigre207 import Cigre207 -from linerate import solver from linerate.equations import ( cigre601, convective_cooling, dimensionless, ieee738, - joule_heating, math, - radiative_cooling, solar_angles, solar_heating, ) from linerate.equations.math import switch_cos_sin +from linerate.models.cigre207 import Cigre207 +from linerate.models.thermal_model import ThermalModel, _copy_method_docstring from linerate.types import Span, Weather from linerate.units import ( Ampere, @@ -39,279 +35,13 @@ __all__ = ["ThermalModel", "Cigre601", "IEEE738", "Cigre207"] -def _copy_method_docstring(parent_class): - def inner(func): - func.__doc__ = getattr(parent_class, func.__name__).__doc__ - return func - - return inner - - -class ThermalModel(ABC): - """Abstract class for a minimal conductor thermal model.""" - - @abstractmethod - def __init__(self, span: Span, weather: Weather): - self.span = span - self.weather = weather - - @abstractmethod - def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: - r"""Compute the conductor resistance, :math:`R~\left[\Omega~\text{m}^{-1}\right]`. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`R~\left[\Omega\right]`. The resistance at the given temperature and current. - """ - resistance = joule_heating.compute_resistance( - conductor_temperature, - temperature1=self.span.conductor.temperature1, - temperature2=self.span.conductor.temperature2, - resistance_at_temperature1=self.span.conductor.resistance_at_temperature1, - resistance_at_temperature2=self.span.conductor.resistance_at_temperature2, - ) - - A = self.span.conductor.aluminium_cross_section_area - b = self.span.conductor.constant_magnetic_effect - m = self.span.conductor.current_density_proportional_magnetic_effect - max_increase = self.span.conductor.max_magnetic_core_relative_resistance_increase - - return joule_heating.correct_resistance_acsr_magnetic_core_loss( - ac_resistance=resistance, - current=current, - aluminium_cross_section_area=A, - constant_magnetic_effect=b, - current_density_proportional_magnetic_effect=m, - max_relative_increase=max_increase, - ) - - @abstractmethod - def compute_joule_heating( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - r"""Compute the Joule heating, :math:`P_J~\left[\text{W}~\text{m}^{-1}\right]`. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_J~\left[\text{W}~\text{m}^{-1}\right]`. The Joule heating. - """ - resistance = self.compute_resistance(conductor_temperature, current) - return joule_heating.compute_joule_heating(current, resistance) - - @abstractmethod - def compute_solar_heating( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - r"""Compute the solar heating, :math:`P_S~\left[\text{W}~\text{m}^{-1}\right]`. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_s~\left[\text{W}~\text{m}^{-1}\right]`. The solar heating. - """ - raise NotImplementedError - - @abstractmethod - def compute_convective_cooling( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - r"""Compute the convective cooling, :math:`P_c~\left[\text{W}~\text{m}^{-1}\right]`. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_c~\left[\text{W}~\text{m}^{-1}\right]`. The convective cooling. - """ - raise NotImplementedError - - @abstractmethod - def compute_radiative_cooling( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - r"""Compute the radiative cooling, :math:`P_r~\left[\text{W}~\text{m}^{-1}\right]`. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_r~\left[\text{W}~\text{m}^{-1}\right]`. The radiative cooling heating. - """ - return radiative_cooling.compute_radiative_cooling( - surface_temperature=conductor_temperature, - air_temperature=self.weather.air_temperature, - conductor_diameter=self.span.conductor.conductor_diameter, - conductor_emissivity=self.span.conductor.emissivity, - ) - - def compute_heat_balance(self, conductor_temperature: Celsius, current: Ampere) -> WattPerMeter: - r"""Compute the conductor's heat balance. Positive means heating and negative means cooling. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`P_J + P_s - P_c - P_r~\left[\text{W}~\text{m}^{-1}\right]`. The heat balance. - """ - P_j = self.compute_joule_heating(conductor_temperature, current) - P_s = self.compute_solar_heating(conductor_temperature, current) - P_c = self.compute_convective_cooling(conductor_temperature, current) - P_r = self.compute_radiative_cooling(conductor_temperature, current) - return P_j + P_s - P_c - P_r - - def compute_info( - self, conductor_temperature: Celsius, current: Ampere - ) -> Dict[str, WattPerMeter]: - r"""Create a dictionary with the different heating and cooling effects. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Dict[str, WattPerMeter] - A dictionary with the magnitude of the different heating and cooling effects. - """ - return { - "convective_cooling": self.compute_convective_cooling(conductor_temperature, current), - "radiative_cooling": self.compute_radiative_cooling(conductor_temperature, current), - "joule_heating": self.compute_joule_heating(conductor_temperature, current), - "solar_heating": self.compute_solar_heating(conductor_temperature, current), - } - - def compute_steady_state_ampacity( - self, - max_conductor_temperature: Celsius, - min_ampacity: Ampere = 0, - max_ampacity: Ampere = 5000, - tolerance: float = 1.0, - ) -> Ampere: - r"""Use the bisection method to compute the steady-state thermal rating (ampacity). - - Parameters - ---------- - max_conductor_temperature: - :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature - min_ampacity: - :math:`I_\text{min}~\left[\text{A}\right]`. Lower bound for the numerical scheme for - computing the ampacity - max_ampacity: - :math:`I_\text{min}~\left[\text{A}\right]`. Upper bound for the numerical scheme for - computing the ampacity - tolerance: - :math:`\Delta I~\left[\text{A}\right]`. The numerical accuracy of the ampacity. The - bisection iterations will stop once the numerical ampacity uncertainty is below - :math:`\Delta I`. The bisection method will run for - :math:`\left\lceil\frac{I_\text{min} - I_\text{min}}{\Delta I}\right\rceil` iterations. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`I~\left[\text{A}\right]`. The thermal rating. - """ - I = solver.compute_conductor_ampacity( # noqa - self.compute_heat_balance, - max_conductor_temperature=max_conductor_temperature, - min_ampacity=min_ampacity, - max_ampacity=max_ampacity, - tolerance=tolerance, - ) - n = self.span.num_conductors - return I * n - - def compute_conductor_temperature( - self, - current: Ampere, - min_temperature: Celsius = -30, - max_temperature: Celsius = 150, - tolerance: float = 0.5, - ) -> Celsius: - r"""Use the bisection method to compute the steady state conductor temperature. - - Parameters - ---------- - current: - :math:`I_\text{max}~\left[\text{A}\right]`. The current flowing through the conductor. - min_temperature: - :math:`T_\text{min}~\left[^\circ\text{C}\right]`. Lower bound for the numerical scheme - for computing the temperature - max_temperature: - :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Upper bound for the numerical scheme - for computing the temperature - tolerance: - :math:`\Delta T~\left[^\circ\text{C}\right]`. The numerical accuracy of the - temperature. The bisection iterations will stop once the numerical temperature - uncertainty is below :math:`\Delta T`. The bisection method will run for - :math:`\left\lceil\frac{T_\text{min} - T_\text{min}}{\Delta T}\right\rceil` iterations. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`I~\left[\text{A}\right]`. The thermal rating. - """ - T = solver.compute_conductor_temperature( - self.compute_heat_balance, - current=current, - min_temperature=min_temperature, - max_temperature=max_temperature, - tolerance=tolerance, - ) - n = self.span.num_conductors - return T / n - - class Cigre601(ThermalModel): def __init__( - self, - span: Span, - weather: Weather, - time: Date, - max_reynolds_number: Real = 4000, # Max value of the angle correction in CIGRE601 + self, + span: Span, + weather: Weather, + time: Date, + max_reynolds_number: Real = 4000, # Max value of the angle correction in CIGRE601 ): super().__init__(span, weather) self.time = time @@ -325,7 +55,7 @@ def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> @_copy_method_docstring(ThermalModel) def compute_joule_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_joule_heating( conductor_temperature=conductor_temperature, current=current @@ -333,7 +63,7 @@ def compute_joule_heating( @_copy_method_docstring(ThermalModel) def compute_solar_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: alpha_s = self.span.conductor.solar_absorptivity F = self.span.ground_albedo @@ -367,7 +97,7 @@ def compute_solar_heating( @_copy_method_docstring(ThermalModel) def compute_convective_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: D = self.span.conductor.conductor_diameter d = self.span.conductor.outer_layer_strand_diameter @@ -423,14 +153,14 @@ def compute_convective_cooling( @_copy_method_docstring(ThermalModel) def compute_radiative_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_radiative_cooling( conductor_temperature=conductor_temperature, current=current ) def compute_temperature_gradient( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> Celsius: r"""Estimate the difference between the core temperature and the surface temperature. @@ -461,11 +191,11 @@ def compute_temperature_gradient( class IEEE738(ThermalModel): def __init__( - self, - span: Span, - weather: Weather, - time: Date, - max_reynolds_number: Real = 50_000, # Max Reynolds number for forced convection + self, + span: Span, + weather: Weather, + time: Date, + max_reynolds_number: Real = 50_000, # Max Reynolds number for forced convection ): super().__init__(span, weather) self.time = time @@ -479,7 +209,7 @@ def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> @_copy_method_docstring(ThermalModel) def compute_joule_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_joule_heating( conductor_temperature=conductor_temperature, current=current @@ -487,7 +217,7 @@ def compute_joule_heating( @_copy_method_docstring(ThermalModel) def compute_solar_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: alpha_s = self.span.conductor.solar_absorptivity # alpha in IEEE phi = self.span.latitude # Lat in IEEE @@ -510,7 +240,7 @@ def compute_solar_heating( @_copy_method_docstring(ThermalModel) def compute_convective_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: D = self.span.conductor.conductor_diameter # D_0 in IEEE y = self.span.conductor_altitude # H_e in IEEE @@ -537,7 +267,7 @@ def compute_convective_cooling( @_copy_method_docstring(ThermalModel) def compute_radiative_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_radiative_cooling( conductor_temperature=conductor_temperature, current=current diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index 96277f2..7d865e6 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -1,6 +1,6 @@ from abc import abstractmethod -from linerate import Span, ThermalModel, Weather +from linerate.types import Span, Weather from linerate.equations import ( cigre207, cigre601, @@ -12,7 +12,7 @@ solar_heating, ) from linerate.equations.math import switch_cos_sin -from linerate.model import _copy_method_docstring +from linerate.models.thermal_model import _copy_method_docstring, ThermalModel from linerate.units import Ampere, Celsius, Date, OhmPerMeter, WattPerMeter diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py new file mode 100644 index 0000000..3ebd139 --- /dev/null +++ b/linerate/models/thermal_model.py @@ -0,0 +1,281 @@ +from abc import ABC, abstractmethod +from typing import Dict + +from linerate import solver +from linerate.equations import ( + joule_heating, + radiative_cooling, +) +from linerate.types import Span, Weather +from linerate.units import ( + Ampere, + Celsius, + OhmPerMeter, + WattPerMeter, +) + + +def _copy_method_docstring(parent_class): + def inner(func): + func.__doc__ = getattr(parent_class, func.__name__).__doc__ + return func + + return inner + + +class ThermalModel(ABC): + """Abstract class for a minimal conductor thermal model.""" + + @abstractmethod + def __init__(self, span: Span, weather: Weather): + self.span = span + self.weather = weather + + @abstractmethod + def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: + r"""Compute the conductor resistance, :math:`R~\left[\Omega~\text{m}^{-1}\right]`. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`R~\left[\Omega\right]`. The resistance at the given temperature and current. + """ + resistance = joule_heating.compute_resistance( + conductor_temperature, + temperature1=self.span.conductor.temperature1, + temperature2=self.span.conductor.temperature2, + resistance_at_temperature1=self.span.conductor.resistance_at_temperature1, + resistance_at_temperature2=self.span.conductor.resistance_at_temperature2, + ) + + A = self.span.conductor.aluminium_cross_section_area + b = self.span.conductor.constant_magnetic_effect + m = self.span.conductor.current_density_proportional_magnetic_effect + max_increase = self.span.conductor.max_magnetic_core_relative_resistance_increase + + return joule_heating.correct_resistance_acsr_magnetic_core_loss( + ac_resistance=resistance, + current=current, + aluminium_cross_section_area=A, + constant_magnetic_effect=b, + current_density_proportional_magnetic_effect=m, + max_relative_increase=max_increase, + ) + + @abstractmethod + def compute_joule_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + r"""Compute the Joule heating, :math:`P_J~\left[\text{W}~\text{m}^{-1}\right]`. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_J~\left[\text{W}~\text{m}^{-1}\right]`. The Joule heating. + """ + resistance = self.compute_resistance(conductor_temperature, current) + return joule_heating.compute_joule_heating(current, resistance) + + @abstractmethod + def compute_solar_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + r"""Compute the solar heating, :math:`P_S~\left[\text{W}~\text{m}^{-1}\right]`. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_s~\left[\text{W}~\text{m}^{-1}\right]`. The solar heating. + """ + raise NotImplementedError + + @abstractmethod + def compute_convective_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + r"""Compute the convective cooling, :math:`P_c~\left[\text{W}~\text{m}^{-1}\right]`. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_c~\left[\text{W}~\text{m}^{-1}\right]`. The convective cooling. + """ + raise NotImplementedError + + @abstractmethod + def compute_radiative_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + r"""Compute the radiative cooling, :math:`P_r~\left[\text{W}~\text{m}^{-1}\right]`. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_r~\left[\text{W}~\text{m}^{-1}\right]`. The radiative cooling heating. + """ + return radiative_cooling.compute_radiative_cooling( + surface_temperature=conductor_temperature, + air_temperature=self.weather.air_temperature, + conductor_diameter=self.span.conductor.conductor_diameter, + conductor_emissivity=self.span.conductor.emissivity, + ) + + def compute_heat_balance(self, conductor_temperature: Celsius, current: Ampere) -> WattPerMeter: + r"""Compute the conductor's heat balance. Positive means heating and negative means cooling. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`P_J + P_s - P_c - P_r~\left[\text{W}~\text{m}^{-1}\right]`. The heat balance. + """ + P_j = self.compute_joule_heating(conductor_temperature, current) + P_s = self.compute_solar_heating(conductor_temperature, current) + P_c = self.compute_convective_cooling(conductor_temperature, current) + P_r = self.compute_radiative_cooling(conductor_temperature, current) + return P_j + P_s - P_c - P_r + + def compute_info( + self, conductor_temperature: Celsius, current: Ampere + ) -> Dict[str, WattPerMeter]: + r"""Create a dictionary with the different heating and cooling effects. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Dict[str, WattPerMeter] + A dictionary with the magnitude of the different heating and cooling effects. + """ + return { + "convective_cooling": self.compute_convective_cooling(conductor_temperature, current), + "radiative_cooling": self.compute_radiative_cooling(conductor_temperature, current), + "joule_heating": self.compute_joule_heating(conductor_temperature, current), + "solar_heating": self.compute_solar_heating(conductor_temperature, current), + } + + def compute_steady_state_ampacity( + self, + max_conductor_temperature: Celsius, + min_ampacity: Ampere = 0, + max_ampacity: Ampere = 5000, + tolerance: float = 1.0, + ) -> Ampere: + r"""Use the bisection method to compute the steady-state thermal rating (ampacity). + + Parameters + ---------- + max_conductor_temperature: + :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature + min_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Lower bound for the numerical scheme for + computing the ampacity + max_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Upper bound for the numerical scheme for + computing the ampacity + tolerance: + :math:`\Delta I~\left[\text{A}\right]`. The numerical accuracy of the ampacity. The + bisection iterations will stop once the numerical ampacity uncertainty is below + :math:`\Delta I`. The bisection method will run for + :math:`\left\lceil\frac{I_\text{min} - I_\text{min}}{\Delta I}\right\rceil` iterations. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I~\left[\text{A}\right]`. The thermal rating. + """ + I = solver.compute_conductor_ampacity( # noqa + self.compute_heat_balance, + max_conductor_temperature=max_conductor_temperature, + min_ampacity=min_ampacity, + max_ampacity=max_ampacity, + tolerance=tolerance, + ) + n = self.span.num_conductors + return I * n + + def compute_conductor_temperature( + self, + current: Ampere, + min_temperature: Celsius = -30, + max_temperature: Celsius = 150, + tolerance: float = 0.5, + ) -> Celsius: + r"""Use the bisection method to compute the steady state conductor temperature. + + Parameters + ---------- + current: + :math:`I_\text{max}~\left[\text{A}\right]`. The current flowing through the conductor. + min_temperature: + :math:`T_\text{min}~\left[^\circ\text{C}\right]`. Lower bound for the numerical scheme + for computing the temperature + max_temperature: + :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Upper bound for the numerical scheme + for computing the temperature + tolerance: + :math:`\Delta T~\left[^\circ\text{C}\right]`. The numerical accuracy of the + temperature. The bisection iterations will stop once the numerical temperature + uncertainty is below :math:`\Delta T`. The bisection method will run for + :math:`\left\lceil\frac{T_\text{min} - T_\text{min}}{\Delta T}\right\rceil` iterations. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I~\left[\text{A}\right]`. The thermal rating. + """ + T = solver.compute_conductor_temperature( + self.compute_heat_balance, + current=current, + min_temperature=min_temperature, + max_temperature=max_temperature, + tolerance=tolerance, + ) + n = self.span.num_conductors + return T / n From 12699b880054346e205d9243bbca5b6fcd19c83c Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Thu, 30 May 2024 13:32:08 +0200 Subject: [PATCH 06/18] Fix formatting --- linerate/model.py | 38 ++++++++++++++++---------------- linerate/models/cigre207.py | 4 ++-- linerate/models/thermal_model.py | 12 ++-------- 3 files changed, 23 insertions(+), 31 deletions(-) diff --git a/linerate/model.py b/linerate/model.py index 9c3aa35..0db8026 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -37,11 +37,11 @@ class Cigre601(ThermalModel): def __init__( - self, - span: Span, - weather: Weather, - time: Date, - max_reynolds_number: Real = 4000, # Max value of the angle correction in CIGRE601 + self, + span: Span, + weather: Weather, + time: Date, + max_reynolds_number: Real = 4000, # Max value of the angle correction in CIGRE601 ): super().__init__(span, weather) self.time = time @@ -55,7 +55,7 @@ def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> @_copy_method_docstring(ThermalModel) def compute_joule_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_joule_heating( conductor_temperature=conductor_temperature, current=current @@ -63,7 +63,7 @@ def compute_joule_heating( @_copy_method_docstring(ThermalModel) def compute_solar_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: alpha_s = self.span.conductor.solar_absorptivity F = self.span.ground_albedo @@ -97,7 +97,7 @@ def compute_solar_heating( @_copy_method_docstring(ThermalModel) def compute_convective_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: D = self.span.conductor.conductor_diameter d = self.span.conductor.outer_layer_strand_diameter @@ -153,14 +153,14 @@ def compute_convective_cooling( @_copy_method_docstring(ThermalModel) def compute_radiative_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_radiative_cooling( conductor_temperature=conductor_temperature, current=current ) def compute_temperature_gradient( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> Celsius: r"""Estimate the difference between the core temperature and the surface temperature. @@ -191,11 +191,11 @@ def compute_temperature_gradient( class IEEE738(ThermalModel): def __init__( - self, - span: Span, - weather: Weather, - time: Date, - max_reynolds_number: Real = 50_000, # Max Reynolds number for forced convection + self, + span: Span, + weather: Weather, + time: Date, + max_reynolds_number: Real = 50_000, # Max Reynolds number for forced convection ): super().__init__(span, weather) self.time = time @@ -209,7 +209,7 @@ def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> @_copy_method_docstring(ThermalModel) def compute_joule_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_joule_heating( conductor_temperature=conductor_temperature, current=current @@ -217,7 +217,7 @@ def compute_joule_heating( @_copy_method_docstring(ThermalModel) def compute_solar_heating( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: alpha_s = self.span.conductor.solar_absorptivity # alpha in IEEE phi = self.span.latitude # Lat in IEEE @@ -240,7 +240,7 @@ def compute_solar_heating( @_copy_method_docstring(ThermalModel) def compute_convective_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: D = self.span.conductor.conductor_diameter # D_0 in IEEE y = self.span.conductor_altitude # H_e in IEEE @@ -267,7 +267,7 @@ def compute_convective_cooling( @_copy_method_docstring(ThermalModel) def compute_radiative_cooling( - self, conductor_temperature: Celsius, current: Ampere + self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: return super().compute_radiative_cooling( conductor_temperature=conductor_temperature, current=current diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index 7d865e6..8af5969 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -1,6 +1,5 @@ from abc import abstractmethod -from linerate.types import Span, Weather from linerate.equations import ( cigre207, cigre601, @@ -12,7 +11,8 @@ solar_heating, ) from linerate.equations.math import switch_cos_sin -from linerate.models.thermal_model import _copy_method_docstring, ThermalModel +from linerate.models.thermal_model import ThermalModel, _copy_method_docstring +from linerate.types import Span, Weather from linerate.units import Ampere, Celsius, Date, OhmPerMeter, WattPerMeter diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py index 3ebd139..fcf004c 100644 --- a/linerate/models/thermal_model.py +++ b/linerate/models/thermal_model.py @@ -2,17 +2,9 @@ from typing import Dict from linerate import solver -from linerate.equations import ( - joule_heating, - radiative_cooling, -) +from linerate.equations import joule_heating, radiative_cooling from linerate.types import Span, Weather -from linerate.units import ( - Ampere, - Celsius, - OhmPerMeter, - WattPerMeter, -) +from linerate.units import Ampere, Celsius, OhmPerMeter, WattPerMeter def _copy_method_docstring(parent_class): From 715b09aa7dacba35e29dfb9a6ec1c0ada6aede98 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 10:42:41 +0200 Subject: [PATCH 07/18] Allow numpy arrays in compute_nusselt_number --- linerate/equations/cigre207/convective_cooling.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py index 3a727d4..ddbc952 100644 --- a/linerate/equations/cigre207/convective_cooling.py +++ b/linerate/equations/cigre207/convective_cooling.py @@ -354,8 +354,7 @@ def compute_nusselt_number( 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 + normal_nusselt = np.maximum(forced_convection_nusselt_number, natural_nusselt_number) + low_wind_nusselt = np.maximum(normal_nusselt, low_wind_nusselt_number) + low_wind_speed = wind_speed < 0.5 + return np.where(low_wind_speed, low_wind_nusselt, normal_nusselt) From 04fa4a2019513d0896dbae57a4bf68346c306a34 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 10:42:58 +0200 Subject: [PATCH 08/18] Minor docs change --- linerate/equations/cigre207/solar_heating.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/linerate/equations/cigre207/solar_heating.py b/linerate/equations/cigre207/solar_heating.py index ef6b4d0..87f1e6c 100644 --- a/linerate/equations/cigre207/solar_heating.py +++ b/linerate/equations/cigre207/solar_heating.py @@ -92,7 +92,7 @@ def compute_global_radiation_intensity( ground). The factor :math:`0.5 \pi` is due the assumption that the ground reflects light diffusely and uniformly in all directions, so the reflected energy is always directed normally to the line. In CIGRE207, it is also assumed that the diffuse radiation is uniformly directed, - which leads to :math:`I_d (0.5 \pi + 0.5 F \pi)` instead of :math:`I_d (1 + 0.5 F \pi)` + which leads to :math:`I_d \pi/2 (1 + F)` instead of :math:`I_d (1 + 0.5 F \pi)` Parameters ---------- From b94b2c9898fb9219084808811dc70c0436bff2de Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 10:43:26 +0200 Subject: [PATCH 09/18] Allow user to disable diffuse radiation --- linerate/models/cigre207.py | 62 +++++++++---------------------------- 1 file changed, 15 insertions(+), 47 deletions(-) diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index 8af5969..87eaf33 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -1,15 +1,10 @@ -from abc import abstractmethod - from linerate.equations import ( cigre207, - cigre601, convective_cooling, dimensionless, - joule_heating, math, solar_angles, - solar_heating, -) + solar_heating, ) from linerate.equations.math import switch_cos_sin from linerate.models.thermal_model import ThermalModel, _copy_method_docstring from linerate.types import Span, Weather @@ -22,9 +17,11 @@ def __init__( span: Span, weather: Weather, time: Date, + include_diffuse_radiation: bool = True ): super().__init__(span, weather) self.time = time + self.include_diffuse_radiation = include_diffuse_radiation @_copy_method_docstring(ThermalModel) def compute_joule_heating( @@ -39,14 +36,12 @@ def compute_solar_heating( self, conductor_temperature: Celsius, current: Ampere ) -> WattPerMeter: alpha_s = self.span.conductor.solar_absorptivity - F = self.span.ground_albedo phi = self.span.latitude gamma_c = self.span.conductor_azimuth y = self.span.conductor_altitude - N_s = self.weather.clearness_ratio D = self.span.conductor.conductor_diameter - omega = solar_angles.compute_hour_angle_relative_to_noon(self.time) + omega = solar_angles.compute_hour_angle_relative_to_noon(self.time, self.span.longitude) delta = solar_angles.compute_solar_declination(self.time) sin_H_s = solar_angles.compute_sin_solar_altitude(phi, delta, omega) chi = solar_angles.compute_solar_azimuth_variable(phi, delta, omega) @@ -57,9 +52,14 @@ def compute_solar_heating( ) sin_eta = switch_cos_sin(cos_eta) - I_B = cigre601.solar_heating.compute_direct_solar_radiation(sin_H_s, N_s, y) - I_d = cigre601.solar_heating.compute_diffuse_sky_radiation(I_B, sin_H_s) - I_T = cigre601.solar_heating.compute_global_radiation_intensity( + I_B = cigre207.solar_heating.compute_direct_solar_radiation(sin_H_s, y) + if self.include_diffuse_radiation: + I_d = cigre207.solar_heating.compute_diffuse_sky_radiation(I_B, sin_H_s) + F = self.span.ground_albedo + else: + I_d = 0 + F = 0 + I_T = cigre207.solar_heating.compute_global_radiation_intensity( I_B, I_d, F, sin_eta, sin_H_s ) return solar_heating.compute_solar_heating( @@ -127,40 +127,8 @@ def compute_radiative_cooling( conductor_temperature=conductor_temperature, current=current ) - @abstractmethod + @_copy_method_docstring(ThermalModel) def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: - r"""Compute the conductor resistance, :math:`R~\left[\Omega~\text{m}^{-1}\right]`. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`R~\left[\Omega\right]`. The resistance at the given temperature and current. - """ - resistance = joule_heating.compute_resistance( - conductor_temperature, - temperature1=self.span.conductor.temperature1, - temperature2=self.span.conductor.temperature2, - resistance_at_temperature1=self.span.conductor.resistance_at_temperature1, - resistance_at_temperature2=self.span.conductor.resistance_at_temperature2, - ) - - A = self.span.conductor.aluminium_cross_section_area - b = self.span.conductor.constant_magnetic_effect - m = self.span.conductor.current_density_proportional_magnetic_effect - max_increase = self.span.conductor.max_magnetic_core_relative_resistance_increase - - return joule_heating.correct_resistance_acsr_magnetic_core_loss( - ac_resistance=resistance, - current=current, - aluminium_cross_section_area=A, - constant_magnetic_effect=b, - current_density_proportional_magnetic_effect=m, - max_relative_increase=max_increase, + return super().compute_resistance( + conductor_temperature=conductor_temperature, current=current ) From 0918b23b7997a89c1e78e4df8f77dad7b268033d Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 12:48:38 +0200 Subject: [PATCH 10/18] Reformat file --- linerate/models/cigre207.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index 87eaf33..4fe5407 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -4,7 +4,8 @@ dimensionless, math, solar_angles, - solar_heating, ) + solar_heating, +) from linerate.equations.math import switch_cos_sin from linerate.models.thermal_model import ThermalModel, _copy_method_docstring from linerate.types import Span, Weather @@ -13,11 +14,7 @@ class Cigre207(ThermalModel): def __init__( - self, - span: Span, - weather: Weather, - time: Date, - include_diffuse_radiation: bool = True + self, span: Span, weather: Weather, time: Date, include_diffuse_radiation: bool = True ): super().__init__(span, weather) self.time = time From 3ae1ae73a11d85b7fe121709a5f95da75137bd30 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 15:05:25 +0200 Subject: [PATCH 11/18] Fix some documentation --- docs/api/equations/math.rst | 2 +- .../equations/ieee738/convective_cooling.py | 8 +- .../test_ratekit_ampacity_cases.py | 115 ++++++++++++++++++ 3 files changed, 120 insertions(+), 5 deletions(-) create mode 100644 tests/acceptance_tests/test_ratekit_ampacity_cases.py diff --git a/docs/api/equations/math.rst b/docs/api/equations/math.rst index ade8725..b615202 100644 --- a/docs/api/equations/math.rst +++ b/docs/api/equations/math.rst @@ -1,5 +1,5 @@ Mathematical utilities ---------------------- -.. automodule:: linerate.equations.cigre601.joule_heating +.. automodule:: linerate.equations.math :members: diff --git a/linerate/equations/ieee738/convective_cooling.py b/linerate/equations/ieee738/convective_cooling.py index aac5e75..58b4ef6 100644 --- a/linerate/equations/ieee738/convective_cooling.py +++ b/linerate/equations/ieee738/convective_cooling.py @@ -27,9 +27,9 @@ def compute_air_temperature_at_boundary_layer( # T_film Parameters ---------- temperature_of_conductor_surface: - :math:'T_s ~left[\circ\text{C}\right]`. The temperature of the surface of the conductor. + :math:`T_s ~left[\circ\text{C}\right]`. The temperature of the surface of the conductor. temperature_of_ambient_air: - :math:'T_a ~left[\circ\text{C}\right]`. The temperature of the ambient air. + :math:`T_a ~left[\circ\text{C}\right]`. The temperature of the ambient air. Returns ------- @@ -275,13 +275,13 @@ def compute_convective_cooling( ) -> WattPerMeter: r"""Compute the convective cooling of the conductor. - On page 11 in :cite:p:`ieee738, it says that one should calculate both forced and natural + On page 11 in :cite:p:`ieee738`, it says that one should calculate both forced and natural convection, and choose the larger of the two as the convective cooling. Parameters ---------- forced_convection: - :math:`q_c. The forced convection. + :math:`q_c`. The forced convection. natural_convection: :math:`q_{cn}~\left[\text{W}~\text{m}^{-1}\right]`. The natural convection. diff --git a/tests/acceptance_tests/test_ratekit_ampacity_cases.py b/tests/acceptance_tests/test_ratekit_ampacity_cases.py new file mode 100644 index 0000000..2e40da6 --- /dev/null +++ b/tests/acceptance_tests/test_ratekit_ampacity_cases.py @@ -0,0 +1,115 @@ +"""Test cases to compare CIGRE TB 207 and IEEE738 to Ratekit.""" + +import numpy as np +import pytest +from pytest import approx + +import linerate + +MILES_PER_M = 1 / 1609.344 +INCH_PER_M = 1 / 0.0254 + + +@pytest.fixture +def parrot_conductor(): + return linerate.Conductor( + core_diameter=0.1004 / INCH_PER_M, + conductor_diameter=1.508 / INCH_PER_M, + outer_layer_strand_diameter=0.1673 / INCH_PER_M, + emissivity=0.8, + solar_absorptivity=0.9, + temperature1=25, + temperature2=75, + resistance_at_temperature2=0.0750 * MILES_PER_M, + resistance_at_temperature1=0.0627 * MILES_PER_M, + aluminium_cross_section_area=float("nan"), + constant_magnetic_effect=1, + current_density_proportional_magnetic_effect=0, + max_magnetic_core_relative_resistance_increase=1, + ) + + +@pytest.fixture +def example_weather_a(): + return linerate.Weather( + air_temperature=np.array(30), + wind_direction=np.radians(0), + wind_speed=0.60, + clearness_ratio=1, + ) + + +@pytest.fixture +def example_span_a(parrot_conductor): + start_tower = linerate.Tower(latitude=65, longitude=0.0000, altitude=0) + end_tower = linerate.Tower(latitude=65, longitude=0.0001, altitude=0) + return linerate.Span( + conductor=parrot_conductor, + start_tower=start_tower, + end_tower=end_tower, + ground_albedo=0.0, + num_conductors=1, + ) + + +@pytest.fixture +def example_model_a(example_span_a, example_weather_a): + return linerate.Cigre207( + example_span_a, + example_weather_a, + np.datetime64("2016-06-21 12:00"), + include_diffuse_radiation=False, + ) + + +def test_example_a_convective_cooling(example_model_a): + assert example_model_a.compute_convective_cooling(50, None) == approx(32.52, abs=0.5) + + +def test_example_a_radiative_cooling(example_model_a): + assert example_model_a.compute_radiative_cooling(50, None) == approx(13.41, abs=0.5) + + +def test_example_a_solar_heating(example_model_a): + assert example_model_a.compute_solar_heating(50, None) == approx(31.04, abs=0.5) + + +def test_example_a_resistance(example_model_a): + assert example_model_a.compute_resistance(50, None) == approx(0.0428 / 1e3, rel=1e-3) + + +def test_example_a_ampacity(example_model_a): + assert example_model_a.compute_steady_state_ampacity( + 50, min_ampacity=0, max_ampacity=10000, tolerance=1e-8 + ) == approx(590, abs=1.5) + + +def test_example_a_convective_cooling_70(example_model_a): + assert example_model_a.compute_convective_cooling(70, None) == approx(63.59, rel=3e-2) + + +def test_example_a_radiative_cooling_70(example_model_a): + assert example_model_a.compute_radiative_cooling(70, None) == approx(29.56, abs=0.5) + + +def test_example_a_solar_heating_70(example_model_a): + assert example_model_a.compute_solar_heating(70, None) == approx(31.05, abs=0.5) + + +def test_example_a_resistance_70(example_model_a): + assert example_model_a.compute_resistance(70, None) == approx(0.0458 / 1e3, rel=1e-3) + + +def test_example_a_ampacity_70(example_model_a): + assert example_model_a.compute_steady_state_ampacity( + 70, min_ampacity=0, max_ampacity=10000, tolerance=1e-8 + ) == approx(1178, abs=1.5) + + +@pytest.fixture +def example_model_b(example_span_a, example_weather_a): + return linerate.IEEE738(example_span_a, example_weather_a, np.datetime64("2016-06-21 12:00")) + + +def test_example_b_solar_heating(example_model_b): + assert example_model_b.compute_solar_heating(50, None) == approx(33.03, abs=0.5) From 3a4e8750cca7283e23b38d9266e9c18617e46e79 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 15:14:58 +0200 Subject: [PATCH 12/18] Add documentation pages for Cigre207 --- docs/api/equations/cigre207/convective_cooling.rst | 5 +++++ docs/api/equations/cigre207/solar_heating.rst | 5 +++++ docs/api/equations/convective_cooling.rst | 5 +++++ docs/api/equations/dimensionless.rst | 5 +++++ docs/api/equations/index.rst | 14 ++++++++++++++ docs/api/equations/solar_heating.rst | 5 +++++ docs/api/model.rst | 6 ++++++ 7 files changed, 45 insertions(+) create mode 100644 docs/api/equations/cigre207/convective_cooling.rst create mode 100644 docs/api/equations/cigre207/solar_heating.rst create mode 100644 docs/api/equations/convective_cooling.rst create mode 100644 docs/api/equations/dimensionless.rst create mode 100644 docs/api/equations/solar_heating.rst diff --git a/docs/api/equations/cigre207/convective_cooling.rst b/docs/api/equations/cigre207/convective_cooling.rst new file mode 100644 index 0000000..127138c --- /dev/null +++ b/docs/api/equations/cigre207/convective_cooling.rst @@ -0,0 +1,5 @@ +Convective cooling in CIGRE207 +------------------------------ + +.. automodule:: linerate.equations.cigre207.convective_cooling + :members: diff --git a/docs/api/equations/cigre207/solar_heating.rst b/docs/api/equations/cigre207/solar_heating.rst new file mode 100644 index 0000000..45dbd17 --- /dev/null +++ b/docs/api/equations/cigre207/solar_heating.rst @@ -0,0 +1,5 @@ +Solar heating in CIGRE207 +------------------------- + +.. automodule:: linerate.equations.cigre207.solar_heating + :members: diff --git a/docs/api/equations/convective_cooling.rst b/docs/api/equations/convective_cooling.rst new file mode 100644 index 0000000..a168efb --- /dev/null +++ b/docs/api/equations/convective_cooling.rst @@ -0,0 +1,5 @@ +Convective cooling +------------------ + +.. automodule:: linerate.equations.convective_cooling + :members: diff --git a/docs/api/equations/dimensionless.rst b/docs/api/equations/dimensionless.rst new file mode 100644 index 0000000..1ca4fc1 --- /dev/null +++ b/docs/api/equations/dimensionless.rst @@ -0,0 +1,5 @@ +Dimensionless numbers +--------------------- + +.. automodule:: linerate.equations.dimensionless + :members: diff --git a/docs/api/equations/index.rst b/docs/api/equations/index.rst index 96b9c99..1a20df9 100644 --- a/docs/api/equations/index.rst +++ b/docs/api/equations/index.rst @@ -13,9 +13,23 @@ Shared functionality :maxdepth: 1 joule_heating + convective_cooling radiative_cooling solar_angles + solar_heating math + dimensionless + +Cigre 207 +^^^^^^^^^ + +.. automodule:: linerate.equations.cigre207 + +.. toctree:: + :maxdepth: 1 + + cigre207/convective_cooling + cigre207/solar_heating Cigre 601 ^^^^^^^^^ diff --git a/docs/api/equations/solar_heating.rst b/docs/api/equations/solar_heating.rst new file mode 100644 index 0000000..9140104 --- /dev/null +++ b/docs/api/equations/solar_heating.rst @@ -0,0 +1,5 @@ +Solar heating +------------- + +.. automodule:: linerate.equations.solar_heating + :members: diff --git a/docs/api/model.rst b/docs/api/model.rst index ca99cd7..ee15270 100644 --- a/docs/api/model.rst +++ b/docs/api/model.rst @@ -8,3 +8,9 @@ The ``model`` module .. autoclass:: linerate.model.Cigre601 :inherited-members: + +.. autoclass:: linerate.model.IEEE738 + :inherited-members: + +.. autoclass:: linerate.model.Cigre207 + :inherited-members: From 24816ffa854bf5d041d6042241b569eedd131bac Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 16:16:42 +0200 Subject: [PATCH 13/18] Fix incorrect symbol --- linerate/equations/cigre207/convective_cooling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py index ddbc952..41e7bcd 100644 --- a/linerate/equations/cigre207/convective_cooling.py +++ b/linerate/equations/cigre207/convective_cooling.py @@ -214,7 +214,7 @@ def compute_low_wind_speed_nusseltnumber( Returns ------- Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`\text{Nu}_{90}`. The corrected Nusselt number for low wind speed. + :math:`\text{Nu}_{cor}`. The corrected Nusselt number for low wind speed. """ return 0.55 * perpendicular_flow_nusselt_number From 8d391b253ce6280ff901da3c6f2f4cbe887be448 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 3 Jun 2024 17:15:45 +0200 Subject: [PATCH 14/18] Remove vectorize and expand natural convection range --- docs/refs.bib | 9 +++ .../equations/cigre207/convective_cooling.py | 77 +++++++------------ 2 files changed, 37 insertions(+), 49 deletions(-) diff --git a/docs/refs.bib b/docs/refs.bib index d7defc1..0b93634 100644 --- a/docs/refs.bib +++ b/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}, diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py index 41e7bcd..3c58ece 100644 --- a/linerate/equations/cigre207/convective_cooling.py +++ b/linerate/equations/cigre207/convective_cooling.py @@ -1,7 +1,6 @@ import warnings import numpy as np -from numba import vectorize from ...units import ( Celsius, @@ -153,7 +152,6 @@ 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, @@ -161,16 +159,11 @@ def _compute_perpendicular_flow_nusseltnumber( # 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 - + 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 @@ -219,24 +212,6 @@ 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, @@ -260,8 +235,17 @@ 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 ) @@ -275,27 +259,24 @@ def _check_horizontal_natural_nusselt_number( 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, ) -> 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( @@ -305,9 +286,7 @@ def compute_horizontal_natural_nusselt_number( 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 ---------- From 5a825bb9d0c50cf3d437bc7d3c04a82f59f80e10 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Tue, 4 Jun 2024 14:38:01 +0200 Subject: [PATCH 15/18] Inline two functions --- .../equations/cigre207/convective_cooling.py | 58 ++++++------------- 1 file changed, 19 insertions(+), 39 deletions(-) diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py index 3c58ece..bad2d82 100644 --- a/linerate/equations/cigre207/convective_cooling.py +++ b/linerate/equations/cigre207/convective_cooling.py @@ -152,21 +152,6 @@ def _check_perpendicular_flow_nusseltnumber_out_of_bounds(reynolds_number): warnings.warn("Reynolds number is out of bounds", stacklevel=5) -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 - 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, @@ -188,10 +173,15 @@ def compute_perpendicular_flow_nusseltnumber( :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, - ) + # From table on page 6 in Cigre207 + Re = reynolds_number + Rs = conductor_roughness + 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 def compute_low_wind_speed_nusseltnumber( @@ -263,22 +253,6 @@ def _check_horizontal_natural_nusselt_number( raise ValueError("GrPr out of bounds: Must be < 10^12.") -def _compute_horizontal_natural_nusselt_number( - grashof_number: Unitless, - prandtl_number: Unitless, -) -> Unitless: - GrPr = grashof_number * prandtl_number - # 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, @@ -301,10 +275,16 @@ def compute_horizontal_natural_nusselt_number( :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, - ) + GrPr = grashof_number * prandtl_number + # 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 def compute_nusselt_number( From 81887bda35774fe0797dba13e8c0706bd1208755 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Tue, 4 Jun 2024 14:47:35 +0200 Subject: [PATCH 16/18] Fix references and symbols --- linerate/equations/cigre207/solar_heating.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/linerate/equations/cigre207/solar_heating.py b/linerate/equations/cigre207/solar_heating.py index 87f1e6c..bcc3750 100644 --- a/linerate/equations/cigre207/solar_heating.py +++ b/linerate/equations/cigre207/solar_heating.py @@ -11,7 +11,7 @@ def compute_direct_solar_radiation( ) -> WattPerSquareMeter: r"""Compute the direct solar radiation. - On page 38 of :cite:p:`cigre601`. Equation (10) states that the direct solar + On page 19 of :cite:p:`cigre601`. Equation (10) states that the direct solar radiation on a surface normal to the solar beam at sea level, :math:`I_{B(0)}`, is given by .. math:: @@ -32,7 +32,7 @@ def compute_direct_solar_radiation( Returns ------- Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. + :math:`I_D~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. """ clearness_ratio = 1.0 return cigre601.solar_heating.compute_direct_solar_radiation( @@ -55,7 +55,7 @@ def compute_diffuse_sky_radiation( Parameters ---------- direct_solar_radiation: - :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. + :math:`I_D~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. sin_solar_altitude: :math:`\sin\left(H_s\right)`. The sine of the solar altitude. @@ -65,8 +65,8 @@ def compute_diffuse_sky_radiation( :math:`I_d~\left[\text{W}~\text{m}^{-2}\right]`.The diffuse solar radiation. """ sin_H_s = sin_solar_altitude - I_B = direct_solar_radiation - return np.maximum(0, (570 - 0.47 * I_B)) * np.maximum(0, sin_H_s) ** 1.2 + I_D = direct_solar_radiation + return np.maximum(0, (570 - 0.47 * I_D)) * np.maximum(0, sin_H_s) ** 1.2 def compute_global_radiation_intensity( @@ -84,7 +84,7 @@ def compute_global_radiation_intensity( .. math:: I_T = - I_B \left(\sin(\eta) + 0.5 F \pi \sin(H_s)\right) + + I_D \left(\sin(\eta) + 0.5 F \pi \sin(H_s)\right) + 0.5 I_d \pi \left(1 + F \right), where :math:`\eta` is the incidence angle of the sun on the line, :math:`H_s` is the solar @@ -97,7 +97,7 @@ def compute_global_radiation_intensity( Parameters ---------- direct_solar_radiation: - :math:`I_B~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. + :math:`I_D~\left[\text{W}~\text{m}^{-2}\right]`. The direct solar radiation. diffuse_sky_radiation: :math:`I_d~\left[\text{W}~\text{m}^{-2}\right]`.The diffuse solar radiation. albedo: @@ -137,11 +137,11 @@ def compute_global_radiation_intensity( * - Snow - 0.6-0.8 """ - I_B = direct_solar_radiation + I_D = direct_solar_radiation I_d = diffuse_sky_radiation F = albedo sin_H_s = sin_solar_altitude sin_eta = sin_angle_of_sun_on_line F_pi_half = 0.5 * pi * F - return I_B * (sin_eta + F_pi_half * sin_H_s) + I_d * pi / 2 * (1 + F) # type: ignore + return I_D * (sin_eta + F_pi_half * sin_H_s) + I_d * pi / 2 * (1 + F) # type: ignore From 0c8afc9e0095a838fafe3a0b44c61c57d340b3b4 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Thu, 6 Jun 2024 17:43:07 +0200 Subject: [PATCH 17/18] Move each model to separate files --- linerate/model.py | 266 +----------------------------------- linerate/models/cigre601.py | 177 ++++++++++++++++++++++++ linerate/models/ieee738.py | 93 +++++++++++++ 3 files changed, 273 insertions(+), 263 deletions(-) create mode 100644 linerate/models/cigre601.py create mode 100644 linerate/models/ieee738.py diff --git a/linerate/model.py b/linerate/model.py index 0db8026..7ba13d1 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -6,269 +6,9 @@ ``linerate.solver`` modules. """ -from numbers import Real - -import numpy as np - -from linerate.equations import ( - cigre601, - convective_cooling, - dimensionless, - ieee738, - math, - solar_angles, - solar_heating, -) -from linerate.equations.math import switch_cos_sin from linerate.models.cigre207 import Cigre207 -from linerate.models.thermal_model import ThermalModel, _copy_method_docstring -from linerate.types import Span, Weather -from linerate.units import ( - Ampere, - Celsius, - Date, - JoulePerKilogramPerKelvin, - OhmPerMeter, - WattPerMeter, -) +from linerate.models.cigre601 import Cigre601 +from linerate.models.ieee738 import IEEE738 +from linerate.models.thermal_model import ThermalModel __all__ = ["ThermalModel", "Cigre601", "IEEE738", "Cigre207"] - - -class Cigre601(ThermalModel): - def __init__( - self, - span: Span, - weather: Weather, - time: Date, - max_reynolds_number: Real = 4000, # Max value of the angle correction in CIGRE601 - ): - super().__init__(span, weather) - self.time = time - self.max_reynolds_number = max_reynolds_number - - @_copy_method_docstring(ThermalModel) - def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: - return super().compute_resistance( - conductor_temperature=conductor_temperature, current=current - ) - - @_copy_method_docstring(ThermalModel) - def compute_joule_heating( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - return super().compute_joule_heating( - conductor_temperature=conductor_temperature, current=current - ) - - @_copy_method_docstring(ThermalModel) - def compute_solar_heating( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - alpha_s = self.span.conductor.solar_absorptivity - F = self.span.ground_albedo - phi = self.span.latitude - gamma_c = self.span.conductor_azimuth - y = self.span.conductor_altitude - N_s = self.weather.clearness_ratio - D = self.span.conductor.conductor_diameter - - omega = solar_angles.compute_hour_angle_relative_to_noon(self.time, self.span.longitude) - delta = solar_angles.compute_solar_declination(self.time) - sin_H_s = solar_angles.compute_sin_solar_altitude(phi, delta, omega) - chi = solar_angles.compute_solar_azimuth_variable(phi, delta, omega) - C = solar_angles.compute_solar_azimuth_constant(chi, omega) - gamma_s = solar_angles.compute_solar_azimuth(C, chi) # Z_c in IEEE - cos_eta = solar_angles.compute_cos_solar_effective_incidence_angle( - sin_H_s, gamma_s, gamma_c - ) - sin_eta = switch_cos_sin(cos_eta) - - I_B = cigre601.solar_heating.compute_direct_solar_radiation(sin_H_s, N_s, y) - I_d = cigre601.solar_heating.compute_diffuse_sky_radiation(I_B, sin_H_s) - I_T = cigre601.solar_heating.compute_global_radiation_intensity( - I_B, I_d, F, sin_eta, sin_H_s - ) - return solar_heating.compute_solar_heating( - alpha_s, - I_T, - D, - ) - - @_copy_method_docstring(ThermalModel) - def compute_convective_cooling( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - D = self.span.conductor.conductor_diameter - d = self.span.conductor.outer_layer_strand_diameter - y = self.span.conductor_altitude - beta = self.span.inclination - V = self.weather.wind_speed - T_a = self.weather.air_temperature - T_c = conductor_temperature - T_f = 0.5 * (T_c + T_a) - - # Compute physical quantities - lambda_f = cigre601.convective_cooling.compute_thermal_conductivity_of_air(T_f) - mu_f = cigre601.convective_cooling.compute_dynamic_viscosity_of_air(T_f) - gamma_f = cigre601.convective_cooling.compute_air_density(T_f, y) - nu_f = cigre601.convective_cooling.compute_kinematic_viscosity_of_air(mu_f, gamma_f) - c_f: JoulePerKilogramPerKelvin = 1005 - delta = math.compute_angle_of_attack( - self.weather.wind_direction, self.span.conductor_azimuth - ) - - # Compute unitless quantities - Re = np.minimum( - dimensionless.compute_reynolds_number(V, D, nu_f), - self.max_reynolds_number, - ) - Gr = dimensionless.compute_grashof_number(D, T_c, T_a, nu_f) - Pr = dimensionless.compute_prandtl_number(lambda_f, mu_f, c_f) - Rs = dimensionless.compute_conductor_roughness(D, d) - - # Compute nusselt numbers - Nu_90 = cigre601.convective_cooling.compute_perpendicular_flow_nusseltnumber( - reynolds_number=Re, conductor_roughness=Rs - ) - Nu_delta = cigre601.convective_cooling.correct_wind_direction_effect_on_nusselt_number( - Nu_90, delta, Rs - ) - - Nu_0 = cigre601.convective_cooling.compute_horizontal_natural_nusselt_number(Gr, Pr) - Nu_beta = cigre601.convective_cooling.correct_natural_nusselt_number_inclination( - Nu_0, beta, Rs - ) - - Nu = cigre601.convective_cooling.compute_nusselt_number( - forced_convection_nusselt_number=Nu_delta, natural_nusselt_number=Nu_beta - ) - - return convective_cooling.compute_convective_cooling( - surface_temperature=conductor_temperature, - air_temperature=self.weather.air_temperature, - nusselt_number=Nu, - thermal_conductivity_of_air=lambda_f, - ) - - @_copy_method_docstring(ThermalModel) - def compute_radiative_cooling( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - return super().compute_radiative_cooling( - conductor_temperature=conductor_temperature, current=current - ) - - def compute_temperature_gradient( - self, conductor_temperature: Celsius, current: Ampere - ) -> Celsius: - r"""Estimate the difference between the core temperature and the surface temperature. - - Parameters - ---------- - conductor_temperature: - :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. - current: - :math:`I~\left[\text{A}\right]`. The current. - - Returns - ------- - Union[float, float64, ndarray[Any, dtype[float64]]] - :math:`T_c - T_s~\left[^\circ \text{C}\right]`. The difference between the core and the - surface temperature of the conductor. - """ - n = self.span.num_conductors - T_c = conductor_temperature - I = current / n # noqa - R = self.compute_resistance(conductor_temperature=T_c, current=I) - return cigre601.convective_cooling.compute_temperature_gradient( - total_heat_gain=I * R, - conductor_thermal_conductivity=self.span.conductor.thermal_conductivity, # type: ignore # noqa - core_diameter=self.span.conductor.core_diameter, - conductor_diameter=self.span.conductor.conductor_diameter, - ) - - -class IEEE738(ThermalModel): - def __init__( - self, - span: Span, - weather: Weather, - time: Date, - max_reynolds_number: Real = 50_000, # Max Reynolds number for forced convection - ): - super().__init__(span, weather) - self.time = time - self.max_reynolds_number = max_reynolds_number - - @_copy_method_docstring(ThermalModel) - def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: - return super().compute_resistance( - conductor_temperature=conductor_temperature, current=current - ) - - @_copy_method_docstring(ThermalModel) - def compute_joule_heating( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - return super().compute_joule_heating( - conductor_temperature=conductor_temperature, current=current - ) - - @_copy_method_docstring(ThermalModel) - def compute_solar_heating( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - alpha_s = self.span.conductor.solar_absorptivity # alpha in IEEE - phi = self.span.latitude # Lat in IEEE - gamma_c = self.span.conductor_azimuth # Z_l i IEEE - y = self.span.conductor_altitude # H_e in IEEE - D = self.span.conductor.conductor_diameter # D_0 in IEEE - - omega = solar_angles.compute_hour_angle_relative_to_noon(self.time, self.span.longitude) - delta = solar_angles.compute_solar_declination(self.time) - sin_H_c = solar_angles.compute_sin_solar_altitude(phi, delta, omega) - Q_s = ieee738.solar_heating.compute_total_heat_flux_density(sin_H_c, True) - K_solar = ieee738.solar_heating.compute_solar_altitude_correction_factor(y) - Q_se = ieee738.solar_heating.compute_elevation_correction_factor(K_solar, Q_s) - chi = solar_angles.compute_solar_azimuth_variable(phi, delta, omega) - C = solar_angles.compute_solar_azimuth_constant(chi, omega) - Z_c = solar_angles.compute_solar_azimuth(C, chi) - cos_theta = solar_angles.compute_cos_solar_effective_incidence_angle(sin_H_c, Z_c, gamma_c) - - return ieee738.solar_heating.compute_solar_heating(alpha_s, Q_se, cos_theta, D) - - @_copy_method_docstring(ThermalModel) - def compute_convective_cooling( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - D = self.span.conductor.conductor_diameter # D_0 in IEEE - y = self.span.conductor_altitude # H_e in IEEE - V = self.weather.wind_speed # V_w in IEEE - T_a = self.weather.air_temperature - T_c = conductor_temperature - T_f = 0.5 * (T_c + T_a) # T_film in IEEE - - mu_f = ieee738.convective_cooling.compute_dynamic_viscosity_of_air(T_f) - rho_f = ieee738.convective_cooling.compute_air_density(T_f, y) - nu_f = ieee738.convective_cooling.compute_kinematic_viscosity_of_air(mu_f, rho_f) - Re = np.minimum( - dimensionless.compute_reynolds_number(V, D, nu_f), # N_Re in IEEE - self.max_reynolds_number, - ) - delta = math.compute_angle_of_attack( - self.weather.wind_direction, self.span.conductor_azimuth - ) # Phi in IEEE - K_angle = ieee738.convective_cooling.compute_wind_direction_factor(delta) - k_f = ieee738.convective_cooling.compute_thermal_conductivity_of_air(T_f) - q_cf = ieee738.convective_cooling.compute_forced_convection(K_angle, Re, k_f, T_c, T_a) - q_cn = ieee738.convective_cooling.compute_natural_convection(rho_f, D, T_c, T_a) - return ieee738.convective_cooling.compute_convective_cooling(q_cf, q_cn) - - @_copy_method_docstring(ThermalModel) - def compute_radiative_cooling( - self, conductor_temperature: Celsius, current: Ampere - ) -> WattPerMeter: - return super().compute_radiative_cooling( - conductor_temperature=conductor_temperature, current=current - ) diff --git a/linerate/models/cigre601.py b/linerate/models/cigre601.py new file mode 100644 index 0000000..5a371ca --- /dev/null +++ b/linerate/models/cigre601.py @@ -0,0 +1,177 @@ +from numbers import Real + +import numpy as np + +from linerate.equations import ( + cigre601, + convective_cooling, + dimensionless, + math, + solar_angles, + solar_heating, +) +from linerate.equations.math import switch_cos_sin +from linerate.models.thermal_model import ThermalModel, _copy_method_docstring +from linerate.types import Span, Weather +from linerate.units import ( + Ampere, + Celsius, + Date, + JoulePerKilogramPerKelvin, + OhmPerMeter, + WattPerMeter, +) + + +class Cigre601(ThermalModel): + def __init__( + self, + span: Span, + weather: Weather, + time: Date, + max_reynolds_number: Real = 4000, # Max value of the angle correction in CIGRE601 + ): + super().__init__(span, weather) + self.time = time + self.max_reynolds_number = max_reynolds_number + + @_copy_method_docstring(ThermalModel) + def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: + return super().compute_resistance( + conductor_temperature=conductor_temperature, current=current + ) + + @_copy_method_docstring(ThermalModel) + def compute_joule_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + return super().compute_joule_heating( + conductor_temperature=conductor_temperature, current=current + ) + + @_copy_method_docstring(ThermalModel) + def compute_solar_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + alpha_s = self.span.conductor.solar_absorptivity + F = self.span.ground_albedo + phi = self.span.latitude + gamma_c = self.span.conductor_azimuth + y = self.span.conductor_altitude + N_s = self.weather.clearness_ratio + D = self.span.conductor.conductor_diameter + + omega = solar_angles.compute_hour_angle_relative_to_noon(self.time, self.span.longitude) + delta = solar_angles.compute_solar_declination(self.time) + sin_H_s = solar_angles.compute_sin_solar_altitude(phi, delta, omega) + chi = solar_angles.compute_solar_azimuth_variable(phi, delta, omega) + C = solar_angles.compute_solar_azimuth_constant(chi, omega) + gamma_s = solar_angles.compute_solar_azimuth(C, chi) # Z_c in IEEE + cos_eta = solar_angles.compute_cos_solar_effective_incidence_angle( + sin_H_s, gamma_s, gamma_c + ) + sin_eta = switch_cos_sin(cos_eta) + + I_B = cigre601.solar_heating.compute_direct_solar_radiation(sin_H_s, N_s, y) + I_d = cigre601.solar_heating.compute_diffuse_sky_radiation(I_B, sin_H_s) + I_T = cigre601.solar_heating.compute_global_radiation_intensity( + I_B, I_d, F, sin_eta, sin_H_s + ) + return solar_heating.compute_solar_heating( + alpha_s, + I_T, + D, + ) + + @_copy_method_docstring(ThermalModel) + def compute_convective_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + D = self.span.conductor.conductor_diameter + d = self.span.conductor.outer_layer_strand_diameter + y = self.span.conductor_altitude + beta = self.span.inclination + V = self.weather.wind_speed + T_a = self.weather.air_temperature + T_c = conductor_temperature + T_f = 0.5 * (T_c + T_a) + + # Compute physical quantities + lambda_f = cigre601.convective_cooling.compute_thermal_conductivity_of_air(T_f) + mu_f = cigre601.convective_cooling.compute_dynamic_viscosity_of_air(T_f) + gamma_f = cigre601.convective_cooling.compute_air_density(T_f, y) + nu_f = cigre601.convective_cooling.compute_kinematic_viscosity_of_air(mu_f, gamma_f) + c_f: JoulePerKilogramPerKelvin = 1005 + delta = math.compute_angle_of_attack( + self.weather.wind_direction, self.span.conductor_azimuth + ) + + # Compute unitless quantities + Re = np.minimum( + dimensionless.compute_reynolds_number(V, D, nu_f), + self.max_reynolds_number, + ) + Gr = dimensionless.compute_grashof_number(D, T_c, T_a, nu_f) + Pr = dimensionless.compute_prandtl_number(lambda_f, mu_f, c_f) + Rs = dimensionless.compute_conductor_roughness(D, d) + + # Compute nusselt numbers + Nu_90 = cigre601.convective_cooling.compute_perpendicular_flow_nusseltnumber( + reynolds_number=Re, conductor_roughness=Rs + ) + Nu_delta = cigre601.convective_cooling.correct_wind_direction_effect_on_nusselt_number( + Nu_90, delta, Rs + ) + + Nu_0 = cigre601.convective_cooling.compute_horizontal_natural_nusselt_number(Gr, Pr) + Nu_beta = cigre601.convective_cooling.correct_natural_nusselt_number_inclination( + Nu_0, beta, Rs + ) + + Nu = cigre601.convective_cooling.compute_nusselt_number( + forced_convection_nusselt_number=Nu_delta, natural_nusselt_number=Nu_beta + ) + + return convective_cooling.compute_convective_cooling( + surface_temperature=conductor_temperature, + air_temperature=self.weather.air_temperature, + nusselt_number=Nu, + thermal_conductivity_of_air=lambda_f, + ) + + @_copy_method_docstring(ThermalModel) + def compute_radiative_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + return super().compute_radiative_cooling( + conductor_temperature=conductor_temperature, current=current + ) + + def compute_temperature_gradient( + self, conductor_temperature: Celsius, current: Ampere + ) -> Celsius: + r"""Estimate the difference between the core temperature and the surface temperature. + + Parameters + ---------- + conductor_temperature: + :math:`T_\text{av}~\left[^\circ\text{C}\right]`. The average conductor temperature. + current: + :math:`I~\left[\text{A}\right]`. The current. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`T_c - T_s~\left[^\circ \text{C}\right]`. The difference between the core and the + surface temperature of the conductor. + """ + n = self.span.num_conductors + T_c = conductor_temperature + I = current / n # noqa + R = self.compute_resistance(conductor_temperature=T_c, current=I) + return cigre601.convective_cooling.compute_temperature_gradient( + total_heat_gain=I * R, + conductor_thermal_conductivity=self.span.conductor.thermal_conductivity, # type: ignore # noqa + core_diameter=self.span.conductor.core_diameter, + conductor_diameter=self.span.conductor.conductor_diameter, + ) diff --git a/linerate/models/ieee738.py b/linerate/models/ieee738.py new file mode 100644 index 0000000..ec5fd98 --- /dev/null +++ b/linerate/models/ieee738.py @@ -0,0 +1,93 @@ +from numbers import Real + +import numpy as np + +from linerate.equations import dimensionless, ieee738, math, solar_angles +from linerate.models.thermal_model import ThermalModel, _copy_method_docstring +from linerate.types import Span, Weather +from linerate.units import Ampere, Celsius, Date, OhmPerMeter, WattPerMeter + + +class IEEE738(ThermalModel): + def __init__( + self, + span: Span, + weather: Weather, + time: Date, + max_reynolds_number: Real = 50_000, # Max Reynolds number for forced convection + ): + super().__init__(span, weather) + self.time = time + self.max_reynolds_number = max_reynolds_number + + @_copy_method_docstring(ThermalModel) + def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: + return super().compute_resistance( + conductor_temperature=conductor_temperature, current=current + ) + + @_copy_method_docstring(ThermalModel) + def compute_joule_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + return super().compute_joule_heating( + conductor_temperature=conductor_temperature, current=current + ) + + @_copy_method_docstring(ThermalModel) + def compute_solar_heating( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + alpha_s = self.span.conductor.solar_absorptivity # alpha in IEEE + phi = self.span.latitude # Lat in IEEE + gamma_c = self.span.conductor_azimuth # Z_l i IEEE + y = self.span.conductor_altitude # H_e in IEEE + D = self.span.conductor.conductor_diameter # D_0 in IEEE + + omega = solar_angles.compute_hour_angle_relative_to_noon(self.time, self.span.longitude) + delta = solar_angles.compute_solar_declination(self.time) + sin_H_c = solar_angles.compute_sin_solar_altitude(phi, delta, omega) + Q_s = ieee738.solar_heating.compute_total_heat_flux_density(sin_H_c, True) + K_solar = ieee738.solar_heating.compute_solar_altitude_correction_factor(y) + Q_se = ieee738.solar_heating.compute_elevation_correction_factor(K_solar, Q_s) + chi = solar_angles.compute_solar_azimuth_variable(phi, delta, omega) + C = solar_angles.compute_solar_azimuth_constant(chi, omega) + Z_c = solar_angles.compute_solar_azimuth(C, chi) + cos_theta = solar_angles.compute_cos_solar_effective_incidence_angle(sin_H_c, Z_c, gamma_c) + + return ieee738.solar_heating.compute_solar_heating(alpha_s, Q_se, cos_theta, D) + + @_copy_method_docstring(ThermalModel) + def compute_convective_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + D = self.span.conductor.conductor_diameter # D_0 in IEEE + y = self.span.conductor_altitude # H_e in IEEE + V = self.weather.wind_speed # V_w in IEEE + T_a = self.weather.air_temperature + T_c = conductor_temperature + T_f = 0.5 * (T_c + T_a) # T_film in IEEE + + mu_f = ieee738.convective_cooling.compute_dynamic_viscosity_of_air(T_f) + rho_f = ieee738.convective_cooling.compute_air_density(T_f, y) + nu_f = ieee738.convective_cooling.compute_kinematic_viscosity_of_air(mu_f, rho_f) + Re = np.minimum( + dimensionless.compute_reynolds_number(V, D, nu_f), # N_Re in IEEE + self.max_reynolds_number, + ) + delta = math.compute_angle_of_attack( + self.weather.wind_direction, self.span.conductor_azimuth + ) # Phi in IEEE + K_angle = ieee738.convective_cooling.compute_wind_direction_factor(delta) + k_f = ieee738.convective_cooling.compute_thermal_conductivity_of_air(T_f) + q_cf = ieee738.convective_cooling.compute_forced_convection(K_angle, Re, k_f, T_c, T_a) + q_cn = ieee738.convective_cooling.compute_natural_convection(rho_f, D, T_c, T_a) + return ieee738.convective_cooling.compute_convective_cooling(q_cf, q_cn) + + @_copy_method_docstring(ThermalModel) + def compute_radiative_cooling( + self, conductor_temperature: Celsius, current: Ampere + ) -> WattPerMeter: + return super().compute_radiative_cooling( + conductor_temperature=conductor_temperature, current=current + ) From 31a81aca193691211ea853c4c647b7e2bdddbd6b Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Mon, 10 Jun 2024 14:54:04 +0200 Subject: [PATCH 18/18] Allow adjustment of direct radiation --- linerate/models/cigre207.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py index 4fe5407..bf95b0d 100644 --- a/linerate/models/cigre207.py +++ b/linerate/models/cigre207.py @@ -14,11 +14,17 @@ class Cigre207(ThermalModel): def __init__( - self, span: Span, weather: Weather, time: Date, include_diffuse_radiation: bool = True + self, + span: Span, + weather: Weather, + time: Date, + include_diffuse_radiation: bool = True, + direct_radiation_factor: float = 1.0, ): super().__init__(span, weather) self.time = time self.include_diffuse_radiation = include_diffuse_radiation + self.direct_radiation_factor = direct_radiation_factor @_copy_method_docstring(ThermalModel) def compute_joule_heating( @@ -49,7 +55,9 @@ def compute_solar_heating( ) sin_eta = switch_cos_sin(cos_eta) - I_B = cigre207.solar_heating.compute_direct_solar_radiation(sin_H_s, y) + I_B = self.direct_radiation_factor * cigre207.solar_heating.compute_direct_solar_radiation( + sin_H_s, y + ) if self.include_diffuse_radiation: I_d = cigre207.solar_heating.compute_diffuse_sky_radiation(I_B, sin_H_s) F = self.span.ground_albedo