Implement temperature tracking and temporary limit calculation

linerate/equations/heat_capacity.py

@@ -0,0 +1,20 @@
+from linerate.units import (
+    Celsius,
+    JoulePerKilogramPerKelvin,
+    KilogramPerMeter,
+    PerKelvin,
+    JoulePerMeter,
+)
+
+
+def calculate_heat_capacity_per_unit_length(
+    mass_per_unit_length: KilogramPerMeter,
+    specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin,
+    specific_heat_capacity_coefficient: PerKelvin,
+    conductor_temperature: Celsius,
+) -> JoulePerMeter:
+    m = mass_per_unit_length
+    c = specific_heat_capacity_at_20_celsius
+    T = conductor_temperature
+    beta = specific_heat_capacity_coefficient
+    return m * c * (1 + (T - 20) * beta)

linerate/models/cigre601.py

@@ -136,9 +136,8 @@ def compute_temperature_gradient(
         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,
+            total_heat_gain=self.compute_joule_heating(conductor_temperature=T_c, current=I),
             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,

linerate/models/thermal_model.py

@@ -1,10 +1,13 @@
 from abc import ABC, abstractmethod
 from typing import Dict
 
+import numpy as np
+from functools import partial
 from linerate import solver
-from linerate.equations import joule_heating, radiative_cooling
-from linerate.types import BaseWeather, Span
-from linerate.units import Ampere, Celsius, OhmPerMeter, WattPerMeter
+from linerate.equations import joule_heating, radiative_cooling, heat_capacity
+from linerate.solver import solve_ivp_forward_euler
+from linerate.types import BaseWeather, ConductorWithTransientData, Span
+from linerate.units import Ampere, Celsius, Duration, OhmPerMeter, WattPerMeter, JoulePerKilogramPerKelvin
 
 
 def _copy_method_docstring(parent_class):
@@ -269,3 +272,126 @@ def compute_conductor_temperature(
             tolerance=tolerance,
         )
         return T
+
+    def compute_heat_capacity_per_unit_length(self, conductor_temperature: Celsius) -> JoulePerKilogramPerKelvin:
+        r"""
+        Compute the heat capacity of the conductor per unit length.
+        Parameters
+        ----------
+        conductor_temperature:
+            :math:`T~\left[^\circ\text{C}\right]` Conductor temperature.
+
+        Returns
+        -------
+        Union[float, float64, ndarray[Any, dtype[float64]]]
+            :math:`c~\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\text{m}^{-1}\right]`. Conductor heat capacity per unit length.
+
+        """
+        conductor = self.span.conductor
+        if not isinstance(conductor, ConductorWithTransientData):
+            raise ValueError("Heat capacity data must be defined to calculate heat capacity.")
+        mc_s = heat_capacity.calculate_heat_capacity_per_unit_length(
+            mass_per_unit_length=conductor.steel_mass_per_unit_length,
+            specific_heat_capacity_at_20_celsius=conductor.steel_specific_heat_capacity_at_20_celsius,
+            specific_heat_capacity_coefficient=conductor.steel_specific_heat_capacity_temperature_coefficient,
+            conductor_temperature=conductor_temperature,
+        )
+        mc_a = heat_capacity.calculate_heat_capacity_per_unit_length(
+            mass_per_unit_length=conductor.aluminum_mass_per_unit_length,
+            specific_heat_capacity_at_20_celsius=conductor.aluminum_specific_heat_capacity_at_20_celsius,
+            specific_heat_capacity_coefficient=conductor.aluminum_specific_heat_capacity_temperature_coefficient,
+            conductor_temperature=conductor_temperature,
+        )
+        return mc_s + mc_a
+
+    def compute_temperature_after_heating(
+        self,
+        initial_conductor_temperature: Celsius,
+        heating_duration: Duration,
+        current: Ampere,
+        time_step: Duration = np.timedelta64(1, "s"),
+    ) -> Celsius:
+        r"""Compute conductor temperature after it is heated for a time using an iterative method.
+
+        Parameters
+        ----------
+        initial_conductor_temperature:
+            :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating.
+        heating_duration:
+            :math:`t`. Duration of heating. Unit is changed to seconds in the function.
+        current:
+            :math:`I~\left[\text{A}\right]`. Current heating the conductor
+        time_step:
+            :math:`\Delta t`. Time step of heating. Unit is changed to seconds in the function.
+
+        Returns
+        -------
+        Union[float, float64, ndarray[Any, dtype[float64]]]
+            :math:`T~\left[^\circ\text{C}\right]`. Conductor temperature.
+        """
+        def conductor_heating(temperature):
+            heat_capacity_ = self.compute_heat_capacity_per_unit_length(temperature)
+            heat_balance = self.compute_heat_balance(temperature, current=current)
+            return heat_balance / heat_capacity_
+        # Convert durations to floating point numbers
+        dt = time_step / np.timedelta64(1, "s")
+        tf = heating_duration / np.timedelta64(1, "s")
+        final_temperature = solve_ivp_forward_euler(conductor_heating, initial_conductor_temperature, tf, dt)
+        return final_temperature
+
+    def compute_transient_ampacity(
+        self,
+        max_conductor_temperature: Celsius,
+        heating_duration: Duration,
+        initial_conductor_temperature: Celsius,
+        time_step: Duration = np.timedelta64(60, "s"),
+        min_ampacity: Ampere = 0,
+        max_ampacity: Ampere = 5000,
+        tolerance: float = 1.0,
+        accept_invalid_values: bool = False,
+    ):
+        r"""Use the bisection method to compute the transient thermal rating (ampacity).
+
+        Parameters
+        ----------
+        max_conductor_temperature:
+            :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature
+        heating_duration:
+            :math:`t`. Duration of heating. Unit is changed to seconds in the function.
+        initial_conductor_temperature:
+            :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating.
+        time_step:
+            :math:`\Delta t`. Time step of heating. Unit is changed to seconds in the function.
+        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.
+        accept_invalid_values:
+            If True, np.nan is returned whenever the current cannot be found within the provided
+            search interval. If False, a ValueError will be raised instead.
+
+
+        Returns
+        -------
+        Union[float, float64, ndarray[Any, dtype[float64]]]
+            :math:`I~\left[\text{A}\right]`. The transient thermal rating.
+        """
+        n = self.span.num_conductors
+        I = solver.compute_conductor_transient_ampacity(  # noqa
+            partial(self.compute_temperature_after_heating, time_step=time_step),
+            max_conductor_temperature,
+            initial_conductor_temperature,
+            heating_duration,
+            min_ampacity,
+            max_ampacity,
+            tolerance,
+            accept_invalid_values,
+        )
+        return I * n

linerate/solver.py

@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from .units import Ampere, Celsius, FloatOrFloatArray, WattPerMeter
+from .units import Ampere, Celsius, Duration, FloatOrFloatArray, WattPerMeter
 
 __all__ = ["bisect", "compute_conductor_temperature", "compute_conductor_ampacity"]
 
@@ -164,3 +164,95 @@ def compute_conductor_ampacity(
     return bisect(
         f, min_ampacity, max_ampacity, tolerance, accept_invalid_values=accept_invalid_values
     )
+
+
+def compute_conductor_transient_ampacity(
+    final_temperature: Callable[[Celsius, Duration, Ampere], WattPerMeter],
+    max_conductor_temperature: Celsius,
+    initial_conductor_temperature: Celsius,
+    heating_duration: Duration,
+    min_ampacity: Ampere = 0,
+    max_ampacity: Ampere = 5_000,
+    tolerance: float = 1,  # Ampere
+    accept_invalid_values: bool = False,
+) -> Ampere:
+    r"""Use the bisection method to compute the temporary thermal rating (ampacity).
+
+    Parameters
+    ----------
+    final_temperature:
+        :math:`f(T, A) = P_J + P_s - P_c - P_r~\left[\text{W}~\text{m}^{-1}\right]`. A function of
+        temperature, heating time and current that returns the conductor temperature at the end of heating time.
+    initial_conductor_temperature:
+        :math:`T_\text{inital}~\left[^\circ\text{C}\right]`. Initial conductor temperature
+    heating_duration:
+        :math:`\delta t_\text{heating}~\left[\text{s}\right]`. Duration the conductor is heated.
+    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{max} - I_\text{min}}{\Delta I}\right\rceil` iterations.
+    accept_invalid_values:
+        If True, np.nan is returned whenever the current cannot be found within the provided
+        search interval. If False, a ValueError will be raised instead.
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+        :math:`I~\left[\text{A}\right]`. The thermal rating.
+    """
+
+    def temperature_difference(current: Ampere) -> Celsius:
+        return max_conductor_temperature - final_temperature(
+            initial_conductor_temperature, heating_duration, current
+        )
+
+    return bisect(
+        temperature_difference,
+        min_ampacity,
+        max_ampacity,
+        tolerance,
+        accept_invalid_values=accept_invalid_values,
+    )
+
+
+def solve_ivp_forward_euler(
+    f: Callable[[FloatOrFloatArray], FloatOrFloatArray],
+    y0: FloatOrFloatArray,
+    duration: FloatOrFloatArray,
+    step: FloatOrFloatArray,
+    ) -> FloatOrFloatArray:
+    """
+    A very basic vectorized forward Euler solver time-invariant systems.
+    Parameters
+    ----------
+    f: Time derivative of the system, assumed to only depend on state. dy / dt = f(y)
+    y0: Initial state of the system
+    duration: Time for which to apply the function
+    step: Time step
+
+    Returns
+    -------
+    Union[float, float64, ndarray[Any, dtype[float64]]]
+    Approximation of function f at the end of the given duration.
+    """
+    y = y0
+    step = np.broadcast_to(step, np.shape(duration))
+    step_count = duration // step
+    remainder = duration % step
+    modification_mask = step_count > 0
+    while np.any(modification_mask):
+        y = y + f(y) * step
+        step_count -= 1
+        modification_mask = step_count > 0
+    if np.any(remainder > 0):
+        y = y + f(y) * remainder
+    return y

linerate/types.py

@@ -8,18 +8,29 @@
 from .units import (
     Celsius,
     Degrees,
+    JoulePerKilogramPerKelvin,
+    KilogramPerMeter,
     Meter,
     MeterPerSecond,
     OhmPerMeter,
+    PerKelvin,
     Radian,
     SquareMeter,
     SquareMeterPerAmpere,
     Unitless,
-    WattPerMeterPerKelvin,
     WattPerSquareMeter,
+    WattPerMeterPerKelvin,
 )
 
-__all__ = ["Conductor", "Weather", "Tower", "Span", "WeatherWithSolarRadiation"]
+__all__ = [
+    "Conductor",
+    "BaseWeather",
+    "Weather",
+    "Tower",
+    "Span",
+    "WeatherWithSolarRadiation",
+    "ConductorWithTransientData",
+]
 
 
 @dataclass(frozen=True)
@@ -75,6 +86,27 @@ class Conductor:
     thermal_conductivity: Optional[WattPerMeterPerKelvin] = None
 
 
+@dataclass(frozen=True, kw_only=True)
+class ConductorWithTransientData(Conductor):
+    #: :math:`m_s`. Mass of steel (core material) per unit length of conductor :math:`\left[\text{kg}\text{m}^{-1}\right]`
+    steel_mass_per_unit_length: KilogramPerMeter
+
+    #: :math:`c_s^20`. Heat capacity of steel (core material) at 20 ^\circ C :math:`\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\right]`
+    steel_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin
+
+    #: :math:`\beta_20^s`. Heat capacity temperature coefficient of the steel (core material) at 20 ^\circ C :math:`\left[text{K}^{-1}\right]`
+    steel_specific_heat_capacity_temperature_coefficient: PerKelvin
+
+    #: :math:`m_a`. Mass of aluminium (outer material) per unit length of conductor :math:`\left[\text{kg}\text{m}^{-1}\right]`
+    aluminum_mass_per_unit_length: KilogramPerMeter
+
+    #: :math:`c_a^20`. Heat capacity of aluminium (outer material) at 20 ^\circ C :math:`\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\right]`
+    aluminum_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin
+
+    #: :math:`\beta_a^20`. Heat capacity temperature coefficient of the aluminium (outer material) at 20 ^\circ C :math:`\left[text{K}^{-1}\right]`
+    aluminum_specific_heat_capacity_temperature_coefficient: PerKelvin
+
+
 @dataclass(frozen=True)
 class Tower:
     """Container for a tower (span end-point)."""

linerate/units.py

@@ -32,5 +32,9 @@
 SquareMeterPerSecond = Annotated[FloatOrFloatArray, "m²/s"]
 KilogramPerMeterPerSecond = Annotated[FloatOrFloatArray, "kg/(m s)"]
 KilogramPerCubeMeter = Annotated[FloatOrFloatArray, "kg/m³"]
+KilogramPerMeter = Annotated[FloatOrFloatArray, "kg/m"]
+PerKelvin = Annotated[FloatOrFloatArray, "1/K"]
+JoulePerMeter = Annotated[FloatOrFloatArray, "J/m"]
 
 Date = Union[np.datetime64, npt.NDArray[np.datetime64]]
+Duration = Union[np.timedelta64, npt.NDArray[np.timedelta64]]

pyproject.toml

@@ -84,7 +84,7 @@ source = ["linerate"]
 log_file = "pytest.log"
 log_level = "DEBUG"
 log_file_format = "%(asctime)s - %(name)s - %(levelname)s: %(message)s"
-norecursedirs = "_build tmp*  __pycache__ src prof wheel_files"
+norecursedirs = ".hypothesis _build tmp*  __pycache__ src prof wheel_files"
 markers = [
   "integration: Mark test as an integration test",
 ]