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/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/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: 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/__init__.py b/linerate/equations/cigre207/__init__.py new file mode 100644 index 0000000..d2c20ed --- /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 diff --git a/linerate/equations/cigre207/ac_resistance.py b/linerate/equations/cigre207/ac_resistance.py new file mode 100644 index 0000000..550128f --- /dev/null +++ b/linerate/equations/cigre207/ac_resistance.py @@ -0,0 +1,25 @@ +from linerate.units import OhmPerMeter + + +def correct_resistance_for_skin_effect( + dc_resistance: OhmPerMeter, +) -> OhmPerMeter: + r""" + Return resistance with constant correction for skin effect, using simple method from + Cigre 207, see section 2.1.1. + + Parameters + ---------- + 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 skin effect into account. + """ + return 1.0123 * dc_resistance + + +# TODO: Implement section 2.1.2? diff --git a/linerate/equations/cigre207/convective_cooling.py b/linerate/equations/cigre207/convective_cooling.py new file mode 100644 index 0000000..bad2d82 --- /dev/null +++ b/linerate/equations/cigre207/convective_cooling.py @@ -0,0 +1,319 @@ +import warnings + +import numpy as np + +from ...units import ( + Celsius, + KilogramPerCubeMeter, + Meter, + MeterPerSecond, + Radian, + SquareMeterPerSecond, + 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 + + +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 +############################# + + +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) + + +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) + # 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( + 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}_{cor}`. The corrected Nusselt number for low wind speed. + """ + return 0.55 * perpendicular_flow_nusselt_number + + +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. + """ + delta = angle_of_attack + Nu_90 = perpendicular_flow_nusselt_number + sin_delta = np.sin(delta) + # Equation (14) on page 7 of Cigre207 + return ( + np.where( + delta <= np.radians(24), + 0.42 + 0.68 * (sin_delta**1.08), + 0.42 + 0.58 * (sin_delta**0.90), + ) + * Nu_90 + ) + + +## Natural convection computations (no wind): +############################################# + + +def _check_horizontal_natural_nusselt_number( + grashof_number: Unitless, prandtl_number: Unitless +) -> None: + GrPr = grashof_number * prandtl_number + if np.any(GrPr < 0): + raise ValueError("GrPr cannot be negative.") + elif np.any(GrPr >= 1e12): + raise ValueError("GrPr out of bounds: Must be < 10^12.") + + +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`. + We expand the allowable range by using Table 9.1 of :cite:p:`incropera2007`. + + 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) + 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( + 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. + """ + 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) diff --git a/linerate/equations/cigre207/solar_heating.py b/linerate/equations/cigre207/solar_heating.py new file mode 100644 index 0000000..bcc3750 --- /dev/null +++ b/linerate/equations/cigre207/solar_heating.py @@ -0,0 +1,147 @@ +import numpy as np +from numpy import pi + +from ...units import Meter, Unitless, WattPerSquareMeter +from .. import cigre601 + + +def compute_direct_solar_radiation( + sin_solar_altitude: Unitless, + height_above_sea_level: Meter, +) -> WattPerSquareMeter: + r"""Compute the direct solar radiation. + + 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:`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. + height_above_sea_level: + :math:`y~\left[\text{m}\right]`. The conductor's altitude. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :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( + sin_solar_altitude, clearness_ratio, height_above_sea_level + ) + + +def compute_diffuse_sky_radiation( + direct_solar_radiation: WattPerSquareMeter, + sin_solar_altitude: Unitless, +) -> WattPerSquareMeter: + r"""Compute the diffuse radiation (light scattered in the atmosphere). + + On page 38 of :cite:p:`cigre207`. + + 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. + + Parameters + ---------- + 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. + + 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_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( + 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 (47) on page 38 of :cite:p:`cigre207` state that the global radiation intensity, + :math:`I_T`, is given by + + .. math:: + + I_T = + 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 + 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 \pi/2 (1 + F)` instead of :math:`I_d (1 + 0.5 F \pi)` + + Parameters + ---------- + 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: + :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_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_D * (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 b4e70bf..2069062 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, @@ -167,169 +163,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 ############################# @@ -634,44 +467,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/cigre601/solar_heating.py b/linerate/equations/cigre601/solar_heating.py index cb2c0bc..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( @@ -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/convective_cooling.py b/linerate/equations/convective_cooling.py new file mode 100644 index 0000000..2c7b427 --- /dev/null +++ b/linerate/equations/convective_cooling.py @@ -0,0 +1,41 @@ +import numpy as np + +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, +) -> 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 diff --git a/linerate/equations/dimensionless.py b/linerate/equations/dimensionless.py new file mode 100644 index 0000000..a1e7c3e --- /dev/null +++ b/linerate/equations/dimensionless.py @@ -0,0 +1,176 @@ +from textwrap import dedent + +import numpy as np + +from linerate.units import ( + Celsius, + JoulePerKilogramPerKelvin, + KilogramPerMeterPerSecond, + Meter, + MeterPerSecond, + MeterPerSquareSecond, + SquareMeterPerSecond, + Unitless, + WattPerMeterPerKelvin, +) + + +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..58b4ef6 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, @@ -28,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 ------- @@ -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: @@ -320,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/linerate/equations/solar_heating.py b/linerate/equations/solar_heating.py new file mode 100644 index 0000000..a177116 --- /dev/null +++ b/linerate/equations/solar_heating.py @@ -0,0 +1,33 @@ +from linerate.units import Meter, Unitless, WattPerMeter, WattPerSquareMeter + + +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 3977c18..7ba13d1 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -6,535 +6,9 @@ ``linerate.solver`` modules. """ -from abc import ABC, abstractmethod -from numbers import Real -from typing import Dict +from linerate.models.cigre207 import Cigre207 +from linerate.models.cigre601 import Cigre601 +from linerate.models.ieee738 import IEEE738 +from linerate.models.thermal_model import ThermalModel -import numpy as np - -from linerate import solver -from linerate.equations import ( - cigre601, - ieee738, - joule_heating, - math, - radiative_cooling, - solar_angles, -) -from linerate.equations.math import switch_cos_sin -from linerate.types import Span, Weather -from linerate.units import ( - Ampere, - Celsius, - Date, - JoulePerKilogramPerKelvin, - OhmPerMeter, - WattPerMeter, -) - -__all__ = ["ThermalModel", "Cigre601", "IEEE738"] - - -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 - ): - 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 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 - 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( - cigre601.convective_cooling.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) - - # 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 cigre601.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( - ieee738.convective_cooling.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 - ) +__all__ = ["ThermalModel", "Cigre601", "IEEE738", "Cigre207"] diff --git a/linerate/models/Cigre207.md b/linerate/models/Cigre207.md new file mode 100644 index 0000000..ad23e73 --- /dev/null +++ b/linerate/models/Cigre207.md @@ -0,0 +1,36 @@ +# 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 temperature 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 for magnetic effects in ACSR +lines as our implementation of Cigre 601. + +## Solar heating + +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, 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. diff --git a/linerate/models/cigre207.py b/linerate/models/cigre207.py new file mode 100644 index 0000000..bf95b0d --- /dev/null +++ b/linerate/models/cigre207.py @@ -0,0 +1,139 @@ +from linerate.equations import ( + cigre207, + 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, OhmPerMeter, WattPerMeter + + +class Cigre207(ThermalModel): + def __init__( + 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( + 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 + phi = self.span.latitude + gamma_c = self.span.conductor_azimuth + y = self.span.conductor_altitude + 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 = 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 + 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( + 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 + 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) + + # 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 + 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) + + # 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 + ) + + @_copy_method_docstring(ThermalModel) + def compute_resistance(self, conductor_temperature: Celsius, current: Ampere) -> OhmPerMeter: + return super().compute_resistance( + 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 + ) diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py new file mode 100644 index 0000000..fcf004c --- /dev/null +++ b/linerate/models/thermal_model.py @@ -0,0 +1,273 @@ +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 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) 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..44695fb --- /dev/null +++ b/tests/equations/cigre207/test_convective_cooling.py @@ -0,0 +1,161 @@ +import numpy as np +from pytest import approx + +from linerate.equations import cigre207, dimensionless + + +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_convective_cooling.py b/tests/equations/cigre601/test_convective_cooling.py index 13a74c4..841f0f8 100644 --- a/tests/equations/cigre601/test_convective_cooling.py +++ b/tests/equations/cigre601/test_convective_cooling.py @@ -259,157 +259,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 +495,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/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/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..b0fc0dc --- /dev/null +++ b/tests/equations/test_convective_cooling.py @@ -0,0 +1,133 @@ +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) + + +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) 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) diff --git a/tests/equations/test_radiative_cooling.py b/tests/equations/test_radiative_cooling.py index 9b8211d..27586bf 100644 --- a/tests/equations/test_radiative_cooling.py +++ b/tests/equations/test_radiative_cooling.py @@ -83,3 +83,16 @@ 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) 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)