"""Attachable subsystems for heat-pump boiler models.
Each subsystem is a self-contained **pure physics engine**: given
physical input state it returns a standardised output dict.
All simulation orchestration (activation logic, result assembly,
exergy calculation) is the responsibility of the scenario class
that uses the subsystem.
Subsystem catalogue
-------------------
- ``SolarThermalCollector`` — flat-plate / evacuated-tube STC
physics engine (``calc_performance``, ``is_preheat_on``)
- ``PhotovoltaicSystem`` — PV + ESS + inverter chain
- ``UVLamp`` — UV disinfection lamp
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import TYPE_CHECKING
import numpy as np
from . import calc_util as cu
from .constants import c_w, k_a, k_D, k_d, rho_w
from .enex_functions import calc_energy_flow
if TYPE_CHECKING:
import pandas as pd
from .dynamic_context import ControlState, StepContext, SubsystemExergy
# ------------------------------------------------------------------
# Default subsystem step() return — STC absent/inactive
# ------------------------------------------------------------------
STC_OFF_STEP: dict = {
"stc_active": False,
"stc_result": {},
"T_stc_w_out_K": np.nan,
"T_stc_pump_w_out_K": np.nan,
"Q_stc_w_out": 0.0,
"Q_stc_pump_w_out": 0.0,
"Q_stc_w_in": 0.0,
"E_stc_pump": 0.0,
"Q_contribution": 0.0,
"E_subsystem": 0.0,
"T_tank_w_in_override_K": None,
}
# Backward-compatible alias
STC_OFF: dict = STC_OFF_STEP
[docs]
@dataclass
class SolarThermalCollector:
"""Solar Thermal Collector (STC) — pure physics engine.
Bundles collector geometry, optical and thermal properties,
and pump parameters. This class is a **stateless physics
calculator**: given physical inputs it returns a performance
dict. All simulation orchestration (activation logic,
result assembly, exergy calculation) is the responsibility
of the scenario class that uses this engine.
The public API consists of:
- :meth:`calc_performance` — single-timestep thermal performance
- :meth:`is_preheat_on` — schedule check
Parameters
----------
A_stc : float
Collector gross area [m²].
stc_tilt : float
Tilt from horizontal [°].
stc_azimuth : float
Azimuth angle (180 = south) [°].
A_stc_pipe : float
Pipe surface area [m²].
alpha_stc : float
Absorptivity [–].
h_o_stc : float
External convective coefficient [W/(m²·K)].
h_r_stc : float
Radiative coefficient [W/(m²·K)].
k_ins_stc : float
Insulation conductivity [W/(m·K)].
x_air_stc : float
Air gap thickness [m].
x_ins_stc : float
Insulation thickness [m].
preheat_start_hour : float
Preheat window start hour.
preheat_end_hour : float
Preheat window end hour.
dV_stc_w : float
Default STC loop flow rate [m³/s].
E_stc_pump : float
STC pump rated power [W].
"""
# Collector
A_stc: float = 2.0
stc_tilt: float = 35.0
stc_azimuth: float = 180.0
A_stc_pipe: float = 2.0
alpha_stc: float = 0.95
h_o_stc: float = 15.0
h_r_stc: float = 2.0
k_ins_stc: float = 0.03
x_air_stc: float = 0.01
x_ins_stc: float = 0.05
# Pump / schedule
preheat_start_hour: float = 6.0
preheat_end_hour: float = 18.0
dV_stc_w: float = 0.001
E_stc_pump: float = 50.0
# ----------------------------------------------------------
# Physics helpers
# ----------------------------------------------------------
[docs]
def calc_overall_heat_transfer_coeff(self) -> float:
"""Compute overall U-value from parallel resistances.
The collector has two heat-loss paths in parallel:
- **Path 1** (top): radiation gap ‖ air gap → external conv
- **Path 2** (bottom): insulation → external conv
Returns
-------
float
Overall heat-loss coefficient [W/(m²·K)].
"""
R_air = self.x_air_stc / k_a
R_ins = self.x_ins_stc / self.k_ins_stc
R_o = 1.0 / self.h_o_stc
R_r = 1.0 / self.h_r_stc
# Path 1: radiation ‖ air gap (parallel), then + external
R1 = (R_r * R_air) / (R_r + R_air) + R_o
# Path 2: insulation (series) + external
R2 = R_ins + R_o
# Two paths in parallel
return 1.0 / R1 + 1.0 / R2
def _calc_solar_absorption(
self,
I_DN_stc: float,
I_dH_stc: float,
) -> tuple[float, float]:
"""Compute total irradiance and absorbed heat.
Parameters
----------
I_DN_stc : float
Direct-normal irradiance [W/m²].
I_dH_stc : float
Diffuse-horizontal irradiance [W/m²].
Returns
-------
tuple[float, float]
``(I_sol_stc, Q_sol_stc)`` — total irradiance and
absorbed solar heat rate [W].
"""
I_sol = I_DN_stc + I_dH_stc
Q_sol = I_sol * self.A_stc_pipe * self.alpha_stc
return I_sol, Q_sol
def _calc_outlet_temperature(
self,
U_stc: float,
G_stc: float,
Q_sol_stc: float,
Q_stc_w_in: float,
T_stc_w_in_K: float,
T0_K: float,
) -> tuple[float, float, float]:
"""Solve for collector outlet and mean plate temps.
Parameters
----------
U_stc : float
Overall U-value [W/(m²·K)].
G_stc : float
Heat capacity flow rate ``c_w · ρ_w · V̇`` [W/K].
Q_sol_stc : float
Absorbed solar heat [W].
Q_stc_w_in : float
Inlet energy flow [W].
T_stc_w_in_K : float
Inlet water temperature [K].
T0_K : float
Dead-state temperature [K].
Returns
-------
tuple[float, float, float]
``(T_stc_w_out_K, T_stc_K, ksi_stc)`` — outlet temp,
mean plate temp, and dimensionless parameter.
"""
A_U = max(1e-6, self.A_stc_pipe * U_stc)
ksi = np.exp(-A_U / G_stc)
T_out_K = ksi * T_stc_w_in_K + (1 - ksi) * T0_K + (1 - ksi) * Q_sol_stc / A_U
T_stc_K = T0_K + (Q_sol_stc - G_stc * (T_out_K - T_stc_w_in_K)) / A_U
return T_out_K, T_stc_K, ksi
# ----------------------------------------------------------
# Main performance calculation
# ----------------------------------------------------------
# ----------------------------------------------------------
# Schedule helper
# ----------------------------------------------------------
[docs]
def is_preheat_on(self, hour_of_day: float) -> bool:
"""Check whether *hour_of_day* falls in the window.
Parameters
----------
hour_of_day : float
Hour within the day (0–24).
Returns
-------
bool
"""
return self.preheat_start_hour <= hour_of_day < self.preheat_end_hour
# ------------------------------------------------------------------
# Photovoltaic System
[docs]
@dataclass
class PhotovoltaicSystem:
"""Photovoltaic System (PV + Charge Controller) — pure physics engine.
Computes PV energy generation from irradiance inputs. This class
is a **stateless physics calculator**: given physical inputs it
returns a performance dict. All routing logic (ESS charge/discharge,
Grid import, dump) is the responsibility of the scenario class.
The public API:
- :meth:`calc_performance` — single-timestep PV generation
"""
# Panel physics
A_pv: float = 5.0
alp_pv: float = 0.9
pv_tilt: float = 35.0
pv_azimuth: float = 180.0
h_o: float = 15.0
# Efficiencies
eta_pv: float = 0.20
eta_ctrl: float = 0.95
T_ctrl_K: float = 308.15
# ------------------------------------------------------------------
# Energy Storage System
# ------------------------------------------------------------------
@dataclass
class EnergyStorageSystem:
"""Energy Storage System (Battery) — pure physics engine.
Accepts charging / discharging *requests* (DC power, [W]) and
enforces capacity and SOC limits. Internal losses and exergy
destruction are computed from round-trip efficiencies.
State (``SOC_ess``) is updated in-place at each call to
:meth:`charge` or :meth:`discharge`. Routing decisions (which
request to fulfil and in what order) are entirely the
responsibility of the scenario class.
Parameters
----------
C_ess_max : float
Rated energy capacity [J]. Default 3.6 MJ (= 1 kWh).
SOC_init : float
Initial state of charge [–].
SOC_min : float
Minimum allowable SOC [–] (depth-of-discharge guard).
SOC_max : float
Maximum allowable SOC [–].
eta_ess_chg : float
Charging efficiency (electricity-in to stored) [–].
eta_ess_dis : float
Discharging efficiency (stored to electricity-out) [–].
T_ess_K : float
Representative battery temperature used for entropy calc [K].
"""
C_ess_max: float = 3.6e6
SOC_init: float = 0.0
SOC_min: float = 0.10
SOC_max: float = 1.00
eta_ess_chg: float = 0.90
eta_ess_dis: float = 0.90
T_ess_K: float = 313.15
SOC_ess: float = field(init=False)
def __post_init__(self) -> None:
self.SOC_ess = self.SOC_init
# ── helpers ──────────────────────────────────────────────────
def _max_chg_power(self, dt: float) -> float:
"""Maximum charging power limited by remaining capacity [W]."""
headroom = (self.SOC_max - self.SOC_ess) * self.C_ess_max
return max(0.0, headroom / (self.eta_ess_chg * dt)) if dt > 0 else 0.0
def _max_dis_power(self, dt: float) -> float:
"""Maximum discharge power limited by SOC_min [W]."""
available = (self.SOC_ess - self.SOC_min) * self.C_ess_max
return max(0.0, available * self.eta_ess_dis / dt) if dt > 0 else 0.0
def _exergy(self, E_chg: float, E_dis: float, T0_K: float) -> dict:
T0 = max(1e-3, T0_K)
Q_l = (1.0 - self.eta_ess_chg) * E_chg + (1.0 / self.eta_ess_dis - 1.0) * E_dis
S_l = Q_l / self.T_ess_K
X_l_ess = Q_l - S_l * T0
X_c_ess = S_l * T0
return {"Q_l_ess": Q_l, "X_l_ess": X_l_ess, "X_c_ess": X_c_ess}
# ── public API ───────────────────────────────────────────────
def charge(self, E_req_chg: float, dt: float, T0_K: float) -> dict:
"""Request to charge the ESS with *E_req_chg* [W] for *dt* [s].
Returns a dict with keys ``E_ess_chg``, ``E_ess_dis``, ``SOC_ess``
plus exergy keys.
"""
E_chg = min(max(E_req_chg, 0.0), self._max_chg_power(dt))
self.SOC_ess += E_chg * self.eta_ess_chg * dt / self.C_ess_max
self.SOC_ess = min(self.SOC_max, self.SOC_ess)
return {
"E_ess_chg": E_chg,
"E_ess_dis": 0.0,
"SOC_ess": self.SOC_ess,
**self._exergy(E_chg, 0.0, T0_K),
}
def discharge(self, E_req_dis: float, dt: float, T0_K: float) -> dict:
"""Request to discharge *E_req_dis* [W] from the ESS for *dt* [s].
Returns a dict with keys ``E_ess_dis`` (actual), ``E_ess_chg``,
``SOC_ess`` plus exergy keys.
"""
E_dis = min(max(E_req_dis, 0.0), self._max_dis_power(dt))
self.SOC_ess -= E_dis / self.eta_ess_dis * dt / self.C_ess_max
self.SOC_ess = max(self.SOC_min, self.SOC_ess)
return {
"E_ess_chg": 0.0,
"E_ess_dis": E_dis,
"SOC_ess": self.SOC_ess,
**self._exergy(0.0, E_dis, T0_K),
}
# ------------------------------------------------------------------
# UV Disinfection Lamp
# ------------------------------------------------------------------
[docs]
@dataclass
class UVLamp:
"""UV disinfection lamp subsystem.
The lamp switches on periodically (``num_switching``
times per ``period_sec``, each for ``exposure_sec``).
All electrical input is converted to heat inside the
tank (``Q_contribution = E_uv``).
Parameters
----------
lamp_watts : float
Rated lamp power [W].
exposure_sec : float
Duration of each on-cycle [s].
num_switching : int
Number of on-cycles per period.
period_sec : float
Switching period [s] (default 3 h = 10 800 s).
"""
lamp_watts: float = 0.0
exposure_sec: float = 0.0
num_switching: int = 1
period_sec: float = 3 * cu.h2s
[docs]
def step(
self,
ctx: StepContext,
ctrl: ControlState,
dt: float,
T_tank_w_in_K: float,
) -> dict:
"""Compute UV lamp state for one timestep."""
from .enex_functions import calc_uv_lamp_power
E_uv: float = calc_uv_lamp_power(
ctx.current_time_s,
self.period_sec,
self.num_switching,
self.exposure_sec,
self.lamp_watts,
)
return {
"Q_contribution": E_uv,
"E_subsystem": E_uv,
"T_tank_w_in_override_K": None,
"E_uv": E_uv,
}
[docs]
def assemble_results(
self,
ctx: StepContext,
ctrl: ControlState,
step_state: dict,
T_solved_K: float,
) -> dict:
"""Report UV power for DataFrame output."""
E_uv: float = step_state.get("E_uv", 0.0)
if E_uv > 0 or self.lamp_watts > 0:
return {"E_uv [W]": E_uv}
return {}
[docs]
def calc_exergy(
self,
df: pd.DataFrame,
T0_K: pd.Series,
) -> SubsystemExergy | None:
"""UV exergy = electricity (handled by E→X conversion).
No additional post-processing needed.
Returns ``None``.
"""
return None