Make simple heat storage works with uniform temperature (FluidMixMapping)

.gitignore

@@ -3,8 +3,12 @@
 build/
 dist/
 .cache/
+.venv
+.env
 .idea
-.cache
+.vscode
+poetry.lock
+.pytest_cache
 *__pycache__*
 *.pyproj
 *.user
@@ -14,4 +18,4 @@ dist/
 doc/.build/*
 doc/_build/*
 doc/build/*
-
+tmp/

doc/source/elements/simple_heat_storage.rst

@@ -43,7 +43,14 @@ Input Static Data
 
     "name", "Custom name for the Storage", "N/A"
     "in_service", "Indicates if the Storage is in service", "N/A"
-    "q_capacity_kwh", "Capacity in kilowatt-hours", "kWh"
+    "e_capacity_kwh", "Capacity in kilowatt-hours (power-only mode)", "kWh"
+    "capacity_kg", "Tank fluid mass; if set, enables FluidMix / uniform tank mode", "kg"
+    "t_tank_init_c", "Initial uniform tank temperature (FluidMix mode)", "°C"
+    "min_temp_c", "Minimum temperature for SOC from T (optional, FluidMix)", "°C"
+    "max_temp_c", "Maximum temperature for SOC from T (optional, FluidMix)", "°C"
+    "u_w_per_m2k", "Wall U-value for heat losses (FluidMix)", "W/(m²·K)"
+    "area_wall_m2", "Wall area for heat losses (FluidMix)", "m²"
+    "t_ext_c", "Ambient temperature for heat losses (FluidMix)", "°C"
 
 
 Input Time Series
@@ -71,7 +78,10 @@ Output Time Series
 Mapping
 ----------
 
-The Simple Storage model model can be mapped using :ref:`GenericMapping <GenericMapping>`.
+The heat storage controller can be connected with:
+
+- **GenericMapping**: power only (input ``q_received_kw``; output ``soc``, ``q_delivered_kw``).
+- **FluidMixMapping**: temperature and mass flows (uniform tank). Set ``capacity_kg`` on the element to enable; optionally set ``t_tank_init_c``, ``min_temp_c``, and ``max_temp_c`` to derive SOC from tank temperature.
 
 
 
@@ -82,7 +92,7 @@ Model
     :members:
 
 
-The heat storage model computes the heat received and delivered by the storage element, and updates the state of charge (SOC) accordingly.
+The heat storage model supports two modes. With **GenericMapping** (power only), it computes the heat received and delivered and updates SOC from the energy balance. With **FluidMixMapping** (element ``capacity_kg`` set), it uses a uniform tank with temperature and mass flows; if ``min_temp_c`` and ``max_temp_c`` are set, SOC is computed from the tank temperature as :math:`\mathrm{SOC} = (T - T_{\min}) / (T_{\max} - T_{\min})` (clipped to [0, 1]).
 
 .. math::
     :nowrap:

pyproject.toml

@@ -33,9 +33,9 @@ dependencies = [
     "converter==1.0.0",
     "matplotlib>=3.9.0",
     "numpy>=1.25.0",
-    "pandapipes>=0.12.0",
+    "pandapipes==0.12.0",
     "pandapower>=2.14.11",
-    "pandas>=1.5.3",
+    "pandas~=2.3",
     "savReaderWriter==3.4.2",
     "scipy>=1.11.0",
     "seaborn==0.12.2",

src/pandaprosumer/constants.py

@@ -2,6 +2,7 @@
 
 CELSIUS_TO_K = NORMAL_TEMPERATURE
 TEMPERATURE_CONVERGENCE_THRESHOLD_C = 1
+MAX_RERUN = 20
 
 
 class HeatExchangerControl:

src/pandaprosumer/controller/base.py

@@ -4,9 +4,11 @@
 
 import logging as pplog
 
+from collections.abc import Iterable
 import numpy as np
 
 from .mapped import MappedController
+from .. import CELSIUS_TO_K
 from ..mapping import FluidMixMapping
 
 logger = pplog.getLogger(__name__)
@@ -48,6 +50,19 @@ def control_step(self, prosumer):
         """
         super().control_step(prosumer)
 
+    def get_cp_fluid_j_per_kgk(self, prosumer, t_c):
+        """
+        Get the heat capacity [J/(kg·K)] of the prosumer's fluid for a temperature t_c [°C].
+        Default to 4180.0 [J/(kg·K)] if no valid fluid is defined in the prosumer.
+        If t_c is a list of temperature, use the average of the temperatures.
+        Use the pandapipes fluid library.
+
+        :param prosumer: The prosumer object
+        :param t_c (float | list[float]): Fluid temperature [°C]
+        :return: float
+        """
+        return prosumer.get_cp_fluid_j_per_kgk(t_c)
+
     def get_treturn_tab_c(self, prosumer):
         """
         Calculates the feed temperature and mass flow to deliver

src/pandaprosumer/controller/data_model/heat_storage.py

@@ -26,4 +26,4 @@ class HeatStorageControllerData:
     element_name: str = 'heat_storage'
     period_index: int = None
     input_columns: List[str] = field(default_factory=lambda: ["q_received_kw"])
-    result_columns: List[str] = field(default_factory=lambda: ["soc", "q_delivered_kw"])
+    result_columns: List[str] = field(default_factory=lambda: ["soc", "t_tank_c", "q_ch_kw", "q_dch_kw", "q_delivered_kw", "mdot_ch_kg_per_s", "t_ch_in_c", "t_ch_out_c", "mdot_dch_kg_per_s", "t_dch_in_c", "t_dch_out_c"])

src/pandaprosumer/controller/models/heat_demand.py

@@ -201,6 +201,7 @@ def control_step(self, prosumer):
             return
 
         if not np.isnan(self._get_input('q_received_kw')):
+            # Generic Mapping case
             q_received_kw = self._get_input('q_received_kw')
             q_uncovered_kw = self._q_demand_kw - q_received_kw
             result = np.array([[q_received_kw, q_uncovered_kw, 0, 0, 0]])
@@ -216,34 +217,45 @@ def control_step(self, prosumer):
             self.applied = True
             return
 
+        # FluidMix Mapping case
         q_demand_kw, t_feed_demand_c, t_return_demand_c, mdot_demand_kg_per_s = self._demand_q_tf_tr_m(prosumer)
 
         # ToDo: If t_in < t_out, return t_in, not t_out
 
-        assert not np.isnan(self._t_in_c), f"Heat Demand {self.name} t_in_c is NaN for timestep {self.time} in prosumer {prosumer.name}"
-        assert not np.isnan(t_feed_demand_c), f"Heat Demand {self.name} t_feed_demand_c is NaN for timestep {self.time} in prosumer {prosumer.name}"
-        assert not np.isnan(t_return_demand_c), f"Heat Demand {self.name} t_return_demand_c is NaN for timestep {self.time} in prosumer {prosumer.name}"
-        assert not np.isnan(q_demand_kw), f"Heat Demand {self.name} q_demand_kw is NaN for timestep {self.time} in prosumer {prosumer.name}"
-        assert not np.isnan(mdot_demand_kg_per_s), f"Heat Demand {self.name} mdot_demand_kg_per_s is NaN for timestep {self.time} in prosumer {prosumer.name}"
-
-        t_mean_c = (self._t_in_c + t_return_demand_c) / 2
-        cp_kj_per_kgk = float(prosumer.fluid.get_heat_capacity(CELSIUS_TO_K + t_mean_c)) / 1000
-        if np.isnan(self._mdot_received_kg_per_s):
+        if np.isnan(self._mdot_received_kg_per_s) and not np.isnan(self._t_in_c):
+            # If upstream sends no flow, but there is a temperature, we can receive some heat with potential 'infinite' masses flow
             q_received_kw = q_demand_kw
+            effective_t_in_c = self._t_in_c
+            t_mean_c = (effective_t_in_c + t_return_demand_c) / 2
+            cp_kj_per_kgk = float(prosumer.fluid.get_heat_capacity(CELSIUS_TO_K + t_mean_c)) / 1000
             mdot_received_kg_per_s = q_demand_kw / (cp_kj_per_kgk * (self._t_in_c - t_return_demand_c))
+            t_out_c = t_return_demand_c
+        elif self._mdot_received_kg_per_s == 0:
+            # No flow from upstream means no heat transfer
+            q_received_kw = 0.0
+            mdot_received_kg_per_s = 0.0
+            effective_t_in_c = self._t_in_c
+            t_out_c = t_return_demand_c  # Should it be effective_t_in_c instead to have t_in_c == t_out_c ?
         else:
+            # Normal operation with flow from upstream
+            effective_t_in_c = self._t_in_c
+            assert not np.isnan(effective_t_in_c), f"Heat Demand {self.name} t_in_c is NaN for timestep {self.time} in prosumer {prosumer.name}"
             mdot_received_kg_per_s = self._mdot_received_kg_per_s
-            q_received_kw = cp_kj_per_kgk * mdot_received_kg_per_s * (self._t_in_c - t_return_demand_c)
-        t_out_c = t_return_demand_c
+
+            t_mean_c = (effective_t_in_c + t_return_demand_c) / 2
+            cp_kj_per_kgk = float(prosumer.fluid.get_heat_capacity(CELSIUS_TO_K + t_mean_c)) / 1000
+            q_received_kw = cp_kj_per_kgk * mdot_received_kg_per_s * (effective_t_in_c - t_return_demand_c)
+            t_out_c = t_return_demand_c
+
         # Calculate the difference between the received and the required power, wo considering the temperature level
         # FixMe: Consider the temperature level in the output
         q_uncovered_kw = q_demand_kw - q_received_kw
-        result = np.array([[q_received_kw, q_uncovered_kw, mdot_received_kg_per_s, self._t_in_c, t_out_c]])
+        result = np.array([[q_received_kw, q_uncovered_kw, mdot_received_kg_per_s, effective_t_in_c, t_out_c]])
         self.last_result = {
             "q_received_kw": q_received_kw,
             "q_uncovered_kw": q_uncovered_kw,
             "mdot_received_kg_per_s": mdot_received_kg_per_s,
-            "t_in_c": self._t_in_c,
+            "t_in_c": effective_t_in_c,
             "t_out_c": t_out_c
         }
         if np.isnan(result).any():

src/pandaprosumer/controller/models/heat_exchanger.py

@@ -7,7 +7,7 @@
 from pandapipes import call_lib
 
 from pandaprosumer.controller.base import BasicProsumerController
-from pandaprosumer.constants import CELSIUS_TO_K, TEMPERATURE_CONVERGENCE_THRESHOLD_C
+from pandaprosumer.constants import CELSIUS_TO_K, TEMPERATURE_CONVERGENCE_THRESHOLD_C, MAX_RERUN
 from pandaprosumer.mapping import FluidMixMapping
 from pandaprosumer.library.heat_exchanger_utils import compute_temp
 
@@ -172,6 +172,13 @@ def calculate_heat_exchanger(self, prosumer, t_2_out_c, t_2_in_c, mdot_2_kg_per_
         delta_t_hot_nom_c = t_1_hot_nom_c - t_2_hot_nom_c
         delta_t_cold_nom_c = t_1_cold_nom_c - t_2_cold_nom_c
         delta_t_hot_c = t_1_in_c - t_2_out_c
+
+        if delta_t_hot_c <= 0:
+            t_1_out_c, mdot_1_kg_per_s = compute_temp(q_ratio, q_exchanged_w, t_1_in_c, t_2_in_c, t_2_out_c,
+                                                      delta_t_hot_nom_c, delta_t_cold_nom_c, cp_1_j_per_kgk,
+                                                      heat_consumer=False)
+            return mdot_1_kg_per_s, t_1_in_c, t_1_out_c, mdot_2_kg_per_s, t_2_in_c, t_2_out_c
+
         # Logarithmic mean temperature difference (LMTD) at nominal conditions
         if delta_t_hot_nom_c == delta_t_cold_nom_c:
             lmtd_nom = delta_t_hot_nom_c
@@ -183,23 +190,42 @@ def calculate_heat_exchanger(self, prosumer, t_2_out_c, t_2_in_c, mdot_2_kg_per_
         max_t_1_out_c = t_1_in_c - min_delta_t_1_c
         min_x = 1 - (max_t_1_out_c - t_2_in_c) / delta_t_hot_c
 
-        min_a = -np.log(1 - min_x) / min_x if min_x != 0 else 1  # FixMe: case where min_x == 1 or >= 1 ?
-
-        if a < min_a:
-            # If 'a' is too low, q_exchanged_w is too big so reduce mdot_2_kg_per_s
-            # else t_1_out_c would be hotter than t_1_in_c
-            a = min_a
-            q_ratio = delta_t_hot_c / (min_a * lmtd_nom)
-            q_exchanged_w = q_ratio * q_exchanged_nom_w
+        if min_x >= 1 :
+            # If t_2_in_c >= max_t_1_out, ignore min_delta_t_1_c but calculate which q_exchanged will give
+            # a higher t_1_out_c
+            # ToDo: create test for this case
+            min_delta_t_cold_c = 3  # ToDo: constant
+            t_1_out_c = t_2_in_c + min_delta_t_cold_c
+            assert t_1_out_c <= t_1_in_c, f"Heat Exchanger {self.name}: t_1_out_c ({t_1_out_c}) is greater than t_1_in_c ({t_1_in_c})"
+            x = 1 - min_delta_t_cold_c / delta_t_hot_c
+            if x > 0 and x < 1:
+                a = -np.log(1 - x) / x
+            else:
+                raise ValueError(f"In prosumer {prosumer.name} a timestep {self.time}: Heat Exchanger {self.name}: Invalid value of x ({x}) calculated from min_delta_t_cold_c ({min_delta_t_cold_c}) and delta_t_hot_c ({delta_t_hot_c})")
+            q_exchanged_w = (delta_t_hot_c * q_exchanged_nom_w) / (a * lmtd_nom)
             mdot_2_kg_per_s = q_exchanged_w / (cp_2_j_per_kgk * delta_t_2_c)
-            # delta_t_cold = _calculate_cold_temperature_difference(a, delta_t_hot_c)
-            t_1_out_c = max_t_1_out_c
             t_mean_1_c = CELSIUS_TO_K + (t_1_in_c + t_1_out_c) / 2
-            mdot_1_kg_per_s = q_exchanged_w / (self.primary_fluid.get_heat_capacity(t_mean_1_c) * (t_1_in_c - t_1_out_c))
+            cp_1_j_per_kgk = self.primary_fluid.get_heat_capacity(t_mean_1_c)
+            mdot_1_kg_per_s = q_exchanged_w / (cp_1_j_per_kgk * (t_1_in_c - t_1_out_c))
+
         else:
-            t_1_out_c, mdot_1_kg_per_s = compute_temp(q_ratio, q_exchanged_w, t_1_in_c, t_2_in_c, t_2_out_c,
-                                                      delta_t_hot_nom_c, delta_t_cold_nom_c, cp_1_j_per_kgk,
-                                                      heat_consumer=False)
+            min_a = -np.log(1 - min_x) / min_x if min_x != 0 else 1
+
+            if a < min_a:
+                # If 'a' is too low, q_exchanged_w is too big so reduce mdot_2_kg_per_s
+                # else t_1_out_c would be hotter than t_1_in_c
+                a = min_a
+                q_ratio = delta_t_hot_c / (min_a * lmtd_nom)
+                q_exchanged_w = q_ratio * q_exchanged_nom_w
+                mdot_2_kg_per_s = q_exchanged_w / (cp_2_j_per_kgk * delta_t_2_c)
+                # delta_t_cold = _calculate_cold_temperature_difference(a, delta_t_hot_c)
+                t_1_out_c = max_t_1_out_c
+                t_mean_1_c = CELSIUS_TO_K + (t_1_in_c + t_1_out_c) / 2
+                mdot_1_kg_per_s = q_exchanged_w / (self.primary_fluid.get_heat_capacity(t_mean_1_c) * (t_1_in_c - t_1_out_c))
+            else:
+                t_1_out_c, mdot_1_kg_per_s = compute_temp(q_ratio, q_exchanged_w, t_1_in_c, t_2_in_c, t_2_out_c,
+                                                          delta_t_hot_nom_c, delta_t_cold_nom_c, cp_1_j_per_kgk,
+                                                          heat_consumer=False)
 
         # If the primary mass flow is too low, no heat is exchanged to the secondary side
         # if mdot_1_kg_per_s < 1e-6:
@@ -297,6 +323,7 @@ def control_step(self, prosumer):
         mdot_2_required_kg_per_s = sum(mdot_tab_required_kg_per_s)
 
         t_1_in_c = self._t_feed_in_c
+
 
         assert not np.isnan(t_1_in_c), f"Heat Exchanger {self.name} t_1_in_c is NaN for timestep {self.time} in prosumer {prosumer.name}"
         assert not np.isnan(t_out_2_required_c), f"Heat Exchanger {self.name} t_out_2_required_c is NaN for timestep {self.time} in prosumer {prosumer.name}"
@@ -312,6 +339,14 @@ def control_step(self, prosumer):
             t_2_out_c = t_out_2_required_c
             result_mdot_tab_kg_per_s = self._merit_order_mass_flow(prosumer, mdot_2_kg_per_s,
                                                                    mdot_tab_required_kg_per_s)
+
+
+            # Handle case where primary mass flow is provided but no heat exchange occurs (full bypass)
+            if not np.isnan(self._mdot_1_provided_kg_per_s):
+                mdot_1_kg_per_s = self._mdot_1_provided_kg_per_s
+                # Temperature remains the same (bypass)
+                t_1_out_c = t_1_in_c
+
         else:
             # FixMe: case where t_out_2_required_c == t_in_2_required_c
             t_out_2_required_c_init = t_out_2_required_c
@@ -327,8 +362,8 @@ def control_step(self, prosumer):
             nb_runs = 0
             while rerun:
                 nb_runs += 1
-                if nb_runs > 20:
-                    raise Exception("Heat Exchanger calculation did not converge after 100 iterations", self.name, self.time, prosumer.name)
+                if nb_runs > MAX_RERUN:
+                    raise Exception(f"Heat Exchanger calculation did not converge after {MAX_RERUN} iterations", self.name, self.time, prosumer.name)
                 (mdot_1_kg_per_s, t_1_in_c, t_1_out_c,
                  mdot_2_kg_per_s, t_2_in_c, t_2_out_c) = self.calculate_heat_exchanger(prosumer,
                                                                                        t_out_2_required_c,

src/pandaprosumer/controller/models/heat_pump.py

@@ -133,7 +133,14 @@ def _calculate_heat_pump(self, prosumer, mdot_cond_kg_per_s, t_cond_out_c, t_con
 
         # 3. Calculate carnot cop
         # FixMe: Why using the condenser output temperature ?
-        cop_carnot = (t_cond_out_c + pinch_c + CELSIUS_TO_K) / (t_cond_out_c - t_evap_in_c)
+        temperature_lift = t_cond_out_c - t_evap_in_c
+        if temperature_lift <= 0:
+            # Heat pump cannot operate if condenser outlet temp is not higher than evaporator inlet temp
+            return (0, 0, 0, 0,
+                    mdot_cond_kg_per_s, t_cond_in_c, t_cond_in_c,
+                    0, t_evap_in_c, t_evap_in_c)
+
+        cop_carnot = (t_cond_out_c + pinch_c + CELSIUS_TO_K) / temperature_lift
 
         # 3bis. Calculate Lorenz cop
         mean_th_c = (t_cond_out_c - t_cond_in_c) / log((t_cond_out_c + CELSIUS_TO_K) / (t_cond_in_c + CELSIUS_TO_K))
@@ -268,7 +275,7 @@ def control_step(self, prosumer):
 
         t_cond_out_required_c, t_cond_in_required_c, mdot_tab_required_kg_per_s = self.t_m_to_deliver(prosumer)
         mdot_cond_required_kg_per_s = np.sum(mdot_tab_required_kg_per_s)
-
+        
         assert not np.isnan(mdot_cond_required_kg_per_s), f"Heat Pump {self.name} mdot_cond_required_kg_per_s is NaN for timestep {self.time} in prosumer {prosumer.name}"
         assert not np.isnan(t_cond_out_required_c), f"Heat Pump {self.name} t_cond_out_required_c is NaN for timestep {self.time} in prosumer {prosumer.name}"
         assert not np.isnan(t_cond_in_required_c), f"Heat Pump {self.name} t_cond_in_required_c is NaN for timestep {self.time} in prosumer {prosumer.name}"