"""Integrated System Model: Air Source Heat Pump Boiler (ASHPB).
This system class orchestrates the dynamic interaction between distinct thermodynamic
sub-components to simulate the overall heating performance. While implemented as an
integrated model, its physical calculations represent the behavior of:
1. **Refrigerant Cycle (Vapor-Compression):**
Evaluates thermodynamic states using CoolProp, enforcing superheat/subcool margins.
2. **Heat Pump Compressor:**
Models the compression process using isentropic and volumetric efficiencies to compute
the actual discharge enthalpy and mass flow rate. The compressor power is determined
from the enthalpy difference and the mass flow rate.
3. **Expansion Valve:**
Modeled as an isenthalpic expansion device (constant enthalpy) that throttles the
refrigerant from the condensing pressure down to the evaporating pressure.
4. **Heat Exchangers (Condenser & Evaporator):**
- **Condenser:** Placed inside the tank (hydronic), utilizing a static overall heat
transfer coefficient (UA_cond_design).
- **Evaporator:** Air-coupled outdoor unit, utilizing a dynamic overall heat transfer
coefficient (UA_evap) that scales non-linearly with fan airflow (Colburn j-factor analogy).
5. **Thermal Storage Tank:**
Modeled with lumped-capacitance and DHW mixing logic.
At each time step, the model finds the minimum-power operating point
(compressor + fan) via bounded 1-D optimisation (Brent's method) over
the evaporator approach temperature difference.
.. note::
이 통합 시스템의 작동 원리와 개별 서브 컴포넌트 모델링에 대한 상세 이론적 배경은
:doc:`/theory/systems/ashp_boiler` 및 :doc:`/theory/components/index` 를 참조하세요.
Optional sub-components (injected via constructor):
- ``SolarThermalCollector`` — tank-circuit or mains-preheat placement
- (future) ``PVPanel`` — photovoltaic integration
Tank-level management and UV disinfection are built-in features
configured through constructor parameters.
"""
import contextlib
import math
from collections.abc import Callable
from dataclasses import dataclass
from typing import Any
import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar
from tqdm import tqdm
from . import calc_util as cu
from .constants import c_a, c_w, rho_a, rho_w
from .dynamic_context import (
ControlState,
StepContext,
determine_heat_source_on_off,
determine_tank_refill_flow,
tank_mass_energy_residual,
)
from .enex_functions import (
calc_HX_perf_for_target_heat,
calc_mixing_valve_flows,
calc_mixing_valve_temp,
)
from .heat_transfer import calc_simple_tank_UA
from .hx_fan import calc_fan_power_from_dV_fan
from .refrigerant import calc_ref_state
from .thermodynamics import calc_energy_flow
from .subsystems import PhotovoltaicSystem, SolarThermalCollector
[docs]
@dataclass
class AirSourceHeatPumpBoiler:
"""Air source heat pump boiler with outdoor-air evaporator.
The refrigerant cycle is resolved via CoolProp with
user-specified superheat / subcool margins. The condenser
approach temperature is determined analytically
(``dT_ref_cond = Q_cond_target / UA_cond``), and a bounded
1-D optimiser (Brent's method) minimises total electrical
input (``E_cmp + E_ou_fan``) over the evaporator approach.
"""
[docs]
def __init__(
self,
# 1. Refrigerant / cycle / compressor -----------
ref: str = "R134a",
V_disp_cmp: float = 0.0002,
eta_cmp_isen: float | Callable | None = None,
eta_cmp_vol: float | Callable | None = None,
eta_cmp_mech: float | Callable = 0.855,
dT_superheat: float = 5.0,
dT_subcool: float = 5.0,
# 2. Heat exchanger -----------------------------
UA_cond_design: float | None = None,
UA_evap_design: float | None = None,
n_evap: float = 0.65,
# 3. Outdoor unit fan ---------------------------
dV_ou_fan_a_design: float | None = None,
dP_ou_fan_design: float = 60.0,
A_cross_ou: float | None = None,
eta_ou_fan_design: float = 0.6,
# 4. Tank / control / load ----------------------
T_tank_w_upper_bound: float = 65.0,
T_tank_w_lower_bound: float = 60.0,
T_mix_w_out: float = 40.0,
T_sup_w: float = 15.0,
hp_capacity: float = 15000.0,
dV_mix_w_out_max: float = 0.0045,
# Tank insulation
r0: float = 0.2,
H: float = 1.2,
x_shell: float = 0.005,
x_ins: float = 0.05,
k_shell: float = 25,
k_ins: float = 0.03,
h_o: float = 15,
# 5. Tank water level management ----------------
tank_always_full: bool = True,
tank_level_lower_bound: float = 0.5,
tank_level_upper_bound: float = 1.0,
dV_tank_w_in_refill: float = 0.001,
prevent_simultaneous_flow: bool = False,
# 7. HP operating schedule ----------------------
hp_on_schedule: list[tuple[float, float]] | None = None,
# 8. Subsystems (class-based injection) ---------
stc: SolarThermalCollector | None = None,
pv: PhotovoltaicSystem | None = None,
uv=None,
# ASHRAE 90.1-2022 VSD coefficients
vsd_coeffs_ou: dict | None = None,
):
if hp_on_schedule is None:
hp_on_schedule = [(0.0, 24.0)]
if vsd_coeffs_ou is None:
vsd_coeffs_ou = {
"c1": 0.0013,
"c2": 0.1470,
"c3": 0.9506,
"c4": -0.0998,
"c5": 0.0,
}
# --- 1. Refrigerant / cycle / compressor ---
self.ref: str = ref
self.V_disp_cmp: float = V_disp_cmp
# Isentropic Efficiency
if eta_cmp_isen is not None:
self.eta_cmp_isen: float | Callable = eta_cmp_isen
else:
self.eta_cmp_isen = 0.80
# Volumetric Efficiency
if eta_cmp_vol is not None:
self.eta_cmp_vol: float | Callable = eta_cmp_vol
else:
self.eta_cmp_vol = lambda r: 0.95 - 0.05 * r
self.eta_cmp_mech: float | Callable = eta_cmp_mech
self.dT_superheat: float = dT_superheat
self.dT_subcool: float = dT_subcool
self.hp_capacity: float = hp_capacity
# --- 2. Heat exchanger UA ---
# If not explicitly provided, the condenser UA dynamically scales to induce a
# ~10.0 K approach temperature difference, which corresponds to the standard
# performance specifications for industrial heat pumps.
# Ref: Application of Industrial Heat Pumps. Annex 35 Final Report (IEA Heat Pump Centre, 2014)
if UA_cond_design is None:
self.UA_cond_design = hp_capacity / 6.0
else:
self.UA_cond_design = UA_cond_design
# The default evaporator UA is determined to ensure an approximate air-side
# temperature drop of 7.0 K across the outdoor unit, aligning with empirical
# laboratory observations of standard residential units.
# Ref: Residential Air Source Heat Pump Water Heater Performance Testing (ORNL, Baxter 2011, DOI: 10.3390/su17052234)
if UA_evap_design is None:
self.UA_evap_design = self.UA_cond_design * 0.8
else:
self.UA_evap_design = UA_evap_design
self.n_evap: float = n_evap
# --- 3. Outdoor unit fan ---
# Default fan flow rate is scaled at 0.0002 m^3/s per W (or 720 CMH per kW),
# representing an optimal ratio of airflow volume to thermal capacity.
# This provides enough margin so that nominal optimization operates at ~80% fan ratio.
if dV_ou_fan_a_design is None:
self.dV_ou_fan_a_design = hp_capacity * 0.00015
else:
self.dV_ou_fan_a_design = dV_ou_fan_a_design
self.dP_ou_fan_design: float = dP_ou_fan_design
self.eta_ou_fan_design: float = eta_ou_fan_design
# Default coil face area assumes a nominal frontal air velocity of 2.0 m/s.
# This velocity is selected specifically to maintain pressure drop profiles
# within optimal ranges for typical plain fin-and-tube configurations.
# Ref: Heat transfer and friction characteristics of plain fin-and-tube heat exchangers, part II (Wang et al., 2000, DOI: 10.1016/S0017-9310(99)00333-6)
if A_cross_ou is None:
self.A_cross_ou = self.dV_ou_fan_a_design / 2.0 # Capped at 2.0 m/s face velocity
else:
self.A_cross_ou = A_cross_ou
self.E_ou_fan_design: float = self.dV_ou_fan_a_design * self.dP_ou_fan_design / self.eta_ou_fan_design
self.vsd_coeffs_ou: dict = vsd_coeffs_ou
self.fan_params_ou: dict = {
"fan_design_flow_rate": self.dV_ou_fan_a_design,
"fan_design_power": self.E_ou_fan_design,
}
# --- 4. Tank geometry and thermal props ---
self.tank_physical: dict = {
"r0": r0,
"H": H,
"x_shell": x_shell,
"x_ins": x_ins,
"k_shell": k_shell,
"k_ins": k_ins,
"h_o": h_o,
}
self.UA_tank: float = calc_simple_tank_UA(
**self.tank_physical,
)
self.V_tank_full: float = math.pi * r0**2 * H
self.C_tank: float = c_w * rho_w * self.V_tank_full
self.dV_mix_w_out_max: float = dV_mix_w_out_max
self.T_tank_w_upper_bound: float = T_tank_w_upper_bound
self.T_tank_w_lower_bound: float = T_tank_w_lower_bound
self.T_sup_w: float = T_sup_w
self.T_sup_w_K: float = cu.C2K(T_sup_w)
self.T_tank_w_in: float = T_sup_w
self.T_mix_w_out: float = T_mix_w_out
self.T_tank_w_in_K: float = cu.C2K(T_sup_w)
self.T_mix_w_out_K: float = cu.C2K(T_mix_w_out)
# --- 5. Tank water level management ---
self.tank_always_full: bool = tank_always_full
self.tank_level_lower_bound: float = tank_level_lower_bound
self.tank_level_upper_bound: float = tank_level_upper_bound
self.dV_tank_w_in_refill: float = dV_tank_w_in_refill
self.prevent_simultaneous_flow: bool = prevent_simultaneous_flow
# --- 6. HP operating schedule ---
self.hp_on_schedule: list[tuple[float, float]] = hp_on_schedule
# --- 7. Subsystems ---
self.stc: SolarThermalCollector | None = stc
self.pv: PhotovoltaicSystem | None = pv
self._subsystems: dict[str, Any] = {}
if stc is not None:
self._subsystems["stc"] = stc
if pv is not None:
self._subsystems["pv"] = pv
if uv is not None:
self._subsystems["uv"] = uv
# Flow-rate sync variables
self.dV_tank_w_in: float = 0.0
self.dV_tank_w_out: float = 0.0
self.dV_mix_sup_w_in: float = 0.0
self.dV_mix_w_out: float = 0.0
# =============================================================
# Refrigerant cycle physics (ASHP-specific)
# =============================================================
def _calc_state(
self,
dT_ref_evap: float,
T_tank_w: float,
Q_ref_cond: float,
T0: float,
*,
flow_state: dict,
) -> dict | None:
"""Evaluate refrigerant cycle at a given operating point.
Parameters
----------
dT_ref_evap : float
Evaporator approach ΔT [K].
T_tank_w : float
Tank water temperature [°C].
Q_ref_cond : float
Target condenser heat rate [W].
T0 : float
Dead-state / outdoor-air temperature [°C].
flow_state : dict
Explicit mixing-valve / tank-flow context. Must contain:
``dV_mix_w_out``, ``dV_tank_w_out``, ``dV_tank_w_in``,
``dV_mix_sup_w_in``. Replaces former implicit ``self.dV_*`` reads.
Returns
-------
dict | None
Cycle performance dictionary; ``None`` if infeasible.
"""
dT_ref_cond: float = Q_ref_cond / self.UA_cond_design if Q_ref_cond > 0 else 0.0
T_tank_w_K: float = cu.C2K(T_tank_w)
T0_K: float = cu.C2K(T0)
T_evap_sat_K: float = T0_K - dT_ref_evap
T_cond_sat_K: float = T_tank_w_K + dT_ref_cond
is_active: bool = Q_ref_cond > 0.0
if not is_active:
# Flow state (explicit parameter, no side-effect reads)
dV_tank_w_out: float = flow_state["dV_tank_w_out [m3/s]"]
dV_tank_w_in: float = flow_state["dV_tank_w_in [m3/s]"]
dV_mix_sup_w_in: float = flow_state["dV_mix_sup_w_in [m3/s]"]
dV_mix_w_out_val: float = flow_state["dV_mix_w_out [m3/s]"]
if dV_mix_w_out_val == 0:
T_mix_w_out_val: float = np.nan
T_mix_w_out_val_K: float = np.nan
else:
mix: dict = calc_mixing_valve_temp(
T_tank_w_K,
self.T_sup_w_K,
self.T_mix_w_out_K,
)
T_mix_w_out_val = mix["T_mix_w_out"]
T_mix_w_out_val_K = mix["T_mix_w_out_K"]
# Energy balance: Q_tank_w_in + Q_ref_cond = Q_tank_w_out + Q_tank_loss + dU_tank/dt
Q_tank_w_in: float = calc_energy_flow(G=c_w * rho_w * dV_tank_w_in, T=self.T_tank_w_in_K, T0=T0_K)
Q_tank_w_out: float = calc_energy_flow(G=c_w * rho_w * dV_tank_w_out, T=T_tank_w_K, T0=T0_K)
Q_mix_sup_w_in: float = calc_energy_flow(G=c_w * rho_w * dV_mix_sup_w_in, T=self.T_sup_w_K, T0=T0_K)
Q_mix_w_out: float = calc_energy_flow(G=c_w * rho_w * dV_mix_w_out_val, T=T_mix_w_out_val_K, T0=T0_K)
cs: dict = calc_ref_state(
T_evap_K=T_evap_sat_K,
T_cond_K=T_cond_sat_K,
refrigerant=self.ref,
eta_cmp_isen=self.eta_cmp_isen,
mode="heating",
dT_superheat=self.dT_superheat,
dT_subcool=0.0,
is_active=False,
)
result: dict = cs.copy()
result.update(
{
"hp_is_on": False,
"converged": True,
# Temperatures [°C]
"T_ou_a_in [°C]": T0,
"T_ou_a_mid [°C]": T0,
"T_ou_a_out [°C]": T0,
"T_tank_w [°C]": T_tank_w,
"T_sup_w [°C]": self.T_sup_w,
"T_tank_w_in [°C]": self.T_tank_w_in,
"T_mix_w_out [°C]": T_mix_w_out_val,
"T0 [°C]": T0,
# Volume flow rates [m3/s]
"dV_ou_a [m3/s]": 0.0,
"v_ou_a [m/s]": 0.0,
"dV_mix_w_out [m3/s]": (dV_mix_w_out_val if dV_mix_w_out_val > 0 else np.nan),
"dV_tank_w_out [m3/s]": (dV_tank_w_out if dV_tank_w_out > 0 else np.nan),
"dV_tank_w_in [m3/s]": (dV_tank_w_in if dV_tank_w_in > 0 else np.nan),
"dV_mix_sup_w_in [m3/s]": (dV_mix_sup_w_in if dV_mix_sup_w_in > 0 else np.nan),
"m_dot_ref [kg/s]": 0.0,
"cmp_rpm [rpm]": 0.0,
# Energy rates [W]
"E_ou_fan [W]": 0.0,
"Q_ref_evap [W]": 0.0,
"Q_ou_a [W]": 0.0,
"E_cmp [W]": 0.0,
"Q_ref_cond [W]": 0.0,
"Q_tank_w_in [W]": Q_tank_w_in,
"Q_tank_w_out [W]": Q_tank_w_out,
"Q_mix_sup_w_in [W]": Q_mix_sup_w_in,
"Q_mix_w_out [W]": Q_mix_w_out,
"E_tot [W]": 0.0,
# COP metrics
"cop_ref [-]": np.nan,
"cop_sys [-]": np.nan,
}
)
return result
# --- Active state calculations ---
pinch_min: float = 0.5
actual_dT_subcool: float = min(self.dT_subcool, max(0.0, dT_ref_cond - pinch_min))
import inspect
def _eval_eff(eff, r_p, rps):
if eff is None:
return 1.0
if callable(eff):
sig = inspect.signature(eff)
if len(sig.parameters) == 2:
return eff(r_p, rps)
return eff(r_p)
return float(eff)
cs: dict = calc_ref_state(
T_evap_K=T_evap_sat_K,
T_cond_K=T_cond_sat_K,
refrigerant=self.ref,
eta_cmp_isen=1.0, # Temporary dummy value to get basic states
mode="heating",
dT_superheat=self.dT_superheat,
dT_subcool=actual_dT_subcool,
is_active=True,
)
ratio_P_cmp = cs["P_ref_cmp_out [Pa]"] / cs["P_ref_cmp_in [Pa]"] if cs["P_ref_cmp_in [Pa]"] > 0 else 1.0
P_cond = cs["P_ref_cmp_out [Pa]"]
s_cmp_in = cs["s_ref_cmp_in [J/(kg·K)]"]
h_cmp_in = cs["h_ref_cmp_in [J/kg]"]
h_exp_in = cs["h_ref_exp_in [J/kg]"]
# Compute isentropic enthalpy once before loop
try:
import CoolProp.CoolProp as CP
h2_isen = CP.PropsSI("H", "P", P_cond, "S", s_cmp_in, self.ref)
except ValueError:
h2_isen = h_cmp_in
def _residual_rps(rps):
val_eta_vol = _eval_eff(self.eta_cmp_vol, ratio_P_cmp, rps)
val_eta_isen = _eval_eff(self.eta_cmp_isen, ratio_P_cmp, rps)
h_cmp_out = h_cmp_in + (h2_isen - h_cmp_in) / val_eta_isen
dh_cond_local = h_cmp_out - h_exp_in
m_dot = self.V_disp_cmp * cs["rho_ref_cmp_in [kg/m3]"] * val_eta_vol * rps
return (m_dot * dh_cond_local) - Q_ref_cond
from scipy.optimize import brentq
try:
cmp_rps = brentq(_residual_rps, 10.0, 150.0)
converged_rps = True
except ValueError:
res_min = _residual_rps(10.0)
res_max = _residual_rps(150.0)
cmp_rps = 10.0 if abs(res_min) < abs(res_max) else 150.0
converged_rps = False
val_eta_vol = _eval_eff(self.eta_cmp_vol, ratio_P_cmp, cmp_rps)
val_eta_isen = _eval_eff(self.eta_cmp_isen, ratio_P_cmp, cmp_rps)
val_eta_mech = _eval_eff(self.eta_cmp_mech, ratio_P_cmp, cmp_rps)
# Post-evaluate state with final isentropic efficiency
cs = calc_ref_state(
T_evap_K=T_evap_sat_K,
T_cond_K=T_cond_sat_K,
refrigerant=self.ref,
eta_cmp_isen=val_eta_isen,
mode="heating",
dT_superheat=self.dT_superheat,
dT_subcool=actual_dT_subcool,
is_active=True,
)
m_dot_ref = self.V_disp_cmp * cs["rho_ref_cmp_in [kg/m3]"] * val_eta_vol * cmp_rps
Q_ref_cond_calc = m_dot_ref * (cs["h_ref_cmp_out [J/kg]"] - cs["h_ref_exp_in [J/kg]"])
Q_ref_evap = m_dot_ref * (cs["h_ref_cmp_in [J/kg]"] - cs["h_ref_exp_out [J/kg]"])
E_cmp = m_dot_ref * (cs["h_ref_cmp_out [J/kg]"] - cs["h_ref_cmp_in [J/kg]"]) / val_eta_mech
HX_perf_ou: dict = calc_HX_perf_for_target_heat(
Q_ref_target=Q_ref_evap,
T_ou_a_in_C=T0,
T_ref_evap_sat_K=cs["T_ref_evap_sat_K"],
T_ref_cond_sat_l_K=cs["T_ref_cond_sat_l_K"],
A_cross=self.A_cross_ou,
UA_design=self.UA_evap_design,
dV_fan_design=self.dV_ou_fan_a_design,
is_active=True,
)
if HX_perf_ou.get("converged", True) is False:
return {
"converged": False,
"_hx_diag": HX_perf_ou,
}
dV_ou_a: float = HX_perf_ou["dV_fan"]
v_ou_a: float = dV_ou_a / self.A_cross_ou
T_ou_a_mid: float = HX_perf_ou["T_ou_a_mid"]
Q_ou_a: float = HX_perf_ou["Q_ou_air"]
E_ou_fan: float = calc_fan_power_from_dV_fan(
dV_fan=dV_ou_a,
fan_params=self.fan_params_ou,
vsd_coeffs=self.vsd_coeffs_ou,
is_active=True,
)
T_ou_a_out: float = T_ou_a_mid + E_ou_fan / (c_a * rho_a * dV_ou_a)
# --- Flow state (explicit parameter, no side-effect reads) ---
dV_tank_w_out: float = flow_state["dV_tank_w_out [m3/s]"]
dV_tank_w_in: float = flow_state["dV_tank_w_in [m3/s]"]
dV_mix_sup_w_in: float = flow_state["dV_mix_sup_w_in [m3/s]"]
dV_mix_w_out_val: float = flow_state["dV_mix_w_out [m3/s]"]
if dV_mix_w_out_val == 0:
T_mix_w_out_val: float = np.nan
T_mix_w_out_val_K: float = np.nan
else:
mix: dict = calc_mixing_valve_temp(
T_tank_w_K,
self.T_sup_w_K,
self.T_mix_w_out_K,
)
T_mix_w_out_val = mix["T_mix_w_out"]
T_mix_w_out_val_K = mix["T_mix_w_out_K"]
# Energy balance: Q_tank_w_in + Q_ref_cond = Q_tank_w_out + Q_tank_loss + dU_tank/dt
Q_tank_w_in: float = calc_energy_flow(G=c_w * rho_w * dV_tank_w_in, T=self.T_tank_w_in_K, T0=T0_K)
Q_tank_w_out: float = calc_energy_flow(G=c_w * rho_w * dV_tank_w_out, T=T_tank_w_K, T0=T0_K)
Q_mix_sup_w_in: float = calc_energy_flow(G=c_w * rho_w * dV_mix_sup_w_in, T=self.T_sup_w_K, T0=T0_K)
Q_mix_w_out: float = calc_energy_flow(G=c_w * rho_w * dV_mix_w_out_val, T=T_mix_w_out_val_K, T0=T0_K)
result: dict = cs.copy()
result.update(
{
"hp_is_on": True,
"converged": bool(converged_rps),
# Temperatures [°C]
"T_ou_a_in [°C]": T0,
"T_ou_a_mid [°C]": T_ou_a_mid,
"T_ou_a_out [°C]": T_ou_a_out,
"T_tank_w [°C]": T_tank_w,
"T_sup_w [°C]": self.T_sup_w,
"T_tank_w_in [°C]": self.T_tank_w_in,
"T_mix_w_out [°C]": T_mix_w_out_val,
"T0 [°C]": T0,
# Volume flow rates [m3/s]
"dV_ou_a [m3/s]": dV_ou_a,
"v_ou_a [m/s]": v_ou_a,
"dV_mix_w_out [m3/s]": (dV_mix_w_out_val if dV_mix_w_out_val > 0 else np.nan),
"dV_tank_w_out [m3/s]": (dV_tank_w_out if dV_tank_w_out > 0 else np.nan),
"dV_tank_w_in [m3/s]": (dV_tank_w_in if dV_tank_w_in > 0 else np.nan),
"dV_mix_sup_w_in [m3/s]": (dV_mix_sup_w_in if dV_mix_sup_w_in > 0 else np.nan),
"m_dot_ref [kg/s]": m_dot_ref, # Mass flow rate [kg/s]
"cmp_rpm [rpm]": cmp_rps * 60, # Compressor speed [rpm]
# Energy rates [W]
"E_ou_fan [W]": E_ou_fan,
"Q_ref_evap [W]": Q_ref_evap,
"Q_ou_a [W]": Q_ou_a,
"E_cmp [W]": E_cmp,
"Q_ref_cond [W]": Q_ref_cond_calc,
"Q_tank_w_in [W]": Q_tank_w_in,
"Q_tank_w_out [W]": Q_tank_w_out,
"Q_mix_sup_w_in [W]": Q_mix_sup_w_in,
"Q_mix_w_out [W]": Q_mix_w_out,
"E_tot [W]": E_cmp + E_ou_fan,
# COP metrics (analogous to X_eff_ref / X_eff_sys)
"cop_ref [-]": (Q_ref_cond_calc / E_cmp if E_cmp > 0 else np.nan),
"cop_sys [-]": (Q_ref_cond_calc / (E_cmp + E_ou_fan) if (E_cmp + E_ou_fan) > 0 else np.nan),
}
)
return result
def _optimize_operation(
self,
T_tank_w: float,
Q_ref_cond: float,
T0: float,
*,
flow_state: dict,
):
"""Find min-power operating point (Brent 1-D).
Parameters
----------
T_tank_w : float
Tank water temperature [°C].
Q_ref_cond : float
Target condenser heat rate [W].
T0 : float
Dead-state temperature [°C].
flow_state : dict
Explicit flow context passed through to ``_calc_state()``.
Returns
-------
scipy.optimize.OptimizeResult
"""
def _objective(dT_ref_evap: float) -> float:
perf: dict | None = self._calc_state(
dT_ref_evap=dT_ref_evap,
T_tank_w=T_tank_w,
Q_ref_cond=Q_ref_cond,
T0=T0,
flow_state=flow_state,
)
if perf is None or not perf.get("converged", False):
return 1e6
E_tot: float = float(perf.get("E_tot [W]", 1e6))
if E_tot <= 0 or np.isnan(E_tot):
return 1e6
return E_tot
return minimize_scalar(
_objective,
bounds=(1.0, 20.0),
method="bounded",
options={"maxiter": 200, "xatol": 1e-6},
)
# =============================================================
# Steady-state analysis
# =============================================================
[docs]
def analyze_steady(
self,
T_tank_w: float,
T0: float,
Q_ref_cond: float,
*,
return_dict: bool = True,
) -> dict | pd.DataFrame:
"""Run a steady-state performance snapshot.
Evaluates the refrigerant cycle at a given operating point
(T_tank_w, T0, Q_ref_cond) **without** solving the tank energy
balance or tracking dynamic flows.
Parameters
----------
T_tank_w : float
Tank water temperature [°C] — treated as a given input.
T0 : float
Dead-state / outdoor-air temperature [°C].
Q_ref_cond : float
Target condenser heat rate [W].
return_dict : bool
If True return dict; else single-row DataFrame.
Returns
-------
dict | pd.DataFrame
"""
import warnings
# Empty flow state as steady state ignores dynamic withdrawal/refill
flow_state = {
"dV_mix_w_out [m3/s]": 0.0,
"dV_tank_w_out [m3/s]": 0.0,
"dV_tank_w_in [m3/s]": 0.0,
"dV_mix_sup_w_in [m3/s]": 0.0,
"alp": 0.0,
}
if Q_ref_cond <= 0:
result = self._calc_state(
dT_ref_evap=5.0,
T_tank_w=T_tank_w,
Q_ref_cond=0.0,
T0=T0,
flow_state=flow_state,
)
else:
opt_result = self._optimize_operation(
T_tank_w=T_tank_w,
Q_ref_cond=Q_ref_cond,
T0=T0,
flow_state=flow_state,
)
result = None
with contextlib.suppress(Exception):
result = self._calc_state(
dT_ref_evap=float(getattr(opt_result, "x", 5.0)),
T_tank_w=T_tank_w,
T0=T0,
Q_ref_cond=Q_ref_cond,
flow_state=flow_state,
)
if result is None or not isinstance(result, dict):
warnings.warn(
f"analyze_steady: optimization failed "
f"(T_tank_w={T_tank_w:.1f}°C, T0={T0:.1f}°C, "
f"Q_ref_cond={Q_ref_cond:.0f}W). "
"Returning HP-off state.",
RuntimeWarning,
stacklevel=2,
)
try:
result = self._calc_state(
dT_ref_evap=5.0,
T_tank_w=T_tank_w,
Q_ref_cond=0.0,
T0=T0,
flow_state=flow_state,
)
except Exception:
result = {
"hp_is_on": False,
"converged": False,
"Q_ref_cond [W]": 0.0,
"Q_ref_evap [W]": 0.0,
"E_cmp [W]": 0.0,
"E_ou_fan [W]": 0.0,
"E_tot [W]": 0.0,
"T_tank_w [°C]": T_tank_w,
"T0 [°C]": T0,
}
if (
result is not None
and isinstance(result, dict)
and "opt_result" in locals()
):
result["converged"] = getattr(opt_result, "success", False)
if result is None or not isinstance(result, dict):
raise RuntimeError("Simulation failed to produce a valid result dictionary.")
# Steady state doesn't have tank loss because we don't solve tank mass/energy balance
result["Q_tank_loss [W]"] = 0.0
result["tank_level [-]"] = 1.0 # steady-state: always_full
if return_dict:
return result
return pd.DataFrame([result])
# =============================================================
# Private helpers for analyze_dynamic
# =============================================================
@staticmethod
def _calc_tank_flow_context(
dV_mix_w_out: float,
T_tank_w_K: float,
T_sup_w_K: float,
T_mix_w_out_K: float,
dV_tank_w_in_override: float | None = None,
) -> dict:
"""Compute mixing-valve / tank-flow context (no side-effects).
Parameters
----------
dV_mix_w_out : float
Service-water draw-off volumetric flow rate [m³/s].
T_tank_w_K : float
Current tank water temperature [K].
T_sup_w_K : float
Mains supply temperature [K].
T_mix_w_out_K : float
Mixing-valve target outlet temperature [K].
dV_tank_w_in_override : float | None
If not None, overrides the symmetry assumption
``dV_tank_w_in = dV_tank_w_out`` (e.g. refill control).
Returns
-------
dict
Keys: ``dV_mix_w_out``, ``dV_tank_w_out``, ``dV_tank_w_in``,
``dV_mix_sup_w_in``.
"""
mix_state = calc_mixing_valve_temp(T_tank_w_K, T_sup_w_K, T_mix_w_out_K)
flows = calc_mixing_valve_flows(dV_mix_w_out, mix_state["alp"])
dV_tank_w_out: float = flows["dV_hot_in"]
dV_tank_w_in: float = dV_tank_w_out if dV_tank_w_in_override is None else dV_tank_w_in_override
return {
"dV_mix_w_out [m3/s]": dV_mix_w_out,
"dV_tank_w_out [m3/s]": dV_tank_w_out,
"dV_tank_w_in [m3/s]": dV_tank_w_in,
"dV_mix_sup_w_in [m3/s]": flows["dV_cold_in"],
"alp": mix_state["alp"]
}
def _determine_hp_state(
self,
ctx: StepContext,
hp_is_on_prev: bool,
) -> tuple[bool, dict, float]:
"""HP on/off + cycle optimisation for one step.
Parameters
----------
ctx : StepContext
Current-step immutable context.
hp_is_on_prev : bool
HP state at previous step.
Returns
-------
tuple[bool, dict, float]
``(hp_is_on, hp_result, Q_ref_cond)``.
"""
T_tank_w: float = cu.K2C(ctx.T_tank_w_K)
hp_is_on: bool = determine_heat_source_on_off(
T_tank_w_C=T_tank_w,
T_lower=self.T_tank_w_lower_bound,
T_upper=self.T_tank_w_upper_bound,
is_on_prev=hp_is_on_prev,
hour_of_day=ctx.hour_of_day,
on_schedule=self.hp_on_schedule,
)
Q_ref_cond: float = self.hp_capacity if hp_is_on else 0.0
# Build explicit flow_state — no side-effects on self.dV_*
flow_state: dict = self._calc_tank_flow_context(
dV_mix_w_out=ctx.dV_mix_w_out,
T_tank_w_K=ctx.T_tank_w_K,
T_sup_w_K=self.T_sup_w_K,
T_mix_w_out_K=self.T_mix_w_out_K,
)
if Q_ref_cond == 0:
hp_result = self._calc_state(
5.0,
T_tank_w,
0.0,
ctx.T0,
flow_state=flow_state,
)
else:
opt = self._optimize_operation(
T_tank_w,
Q_ref_cond,
ctx.T0,
flow_state=flow_state,
)
hp_result = self._calc_state(
float(getattr(opt, "x", 5.0)),
T_tank_w,
Q_ref_cond,
ctx.T0,
flow_state=flow_state,
)
if hp_result is None:
hp_result = {}
return (
hp_is_on,
hp_result,
float(hp_result.get("Q_ref_cond [W]", 0.0)),
)
def _assemble_core_results(
self,
ctx: StepContext,
ctrl: ControlState,
T_solved_K: float,
level_solved: float,
ier: int,
) -> dict:
"""Build HP-core result dict at solved state.
Subsystem results are appended separately by
each subsystem's ``assemble_results()``.
"""
den: float = max(
1e-6,
T_solved_K - self.T_sup_w_K,
)
alp: float = min(
1.0,
max(
0.0,
(self.T_mix_w_out_K - self.T_sup_w_K) / den,
),
)
dV_tank_w_out: float = alp * ctx.dV_mix_w_out
dV_tank_w_in: float = dV_tank_w_out if ctrl.dV_tank_w_in_ctrl is None else ctrl.dV_tank_w_in_ctrl
self.dV_tank_w_out = dV_tank_w_out
self.dV_tank_w_in = dV_tank_w_in
self.dV_mix_w_out = ctx.dV_mix_w_out
self.dV_mix_sup_w_in = (1 - alp) * ctx.dV_mix_w_out
T_mix_w_out_val: float = (
calc_mixing_valve_temp(
T_solved_K,
self.T_sup_w_K,
self.T_mix_w_out_K,
)["T_mix_w_out"]
if ctx.dV_mix_w_out > 0
else np.nan
)
r: dict = {}
r.update(ctrl.result)
r.update(
{
"hp_is_on": ctrl.is_on,
"Q_tank_loss [W]": (self.UA_tank * (T_solved_K - ctx.T0_K)),
"T_tank_w [°C]": cu.K2C(T_solved_K),
"T_mix_w_out [°C]": T_mix_w_out_val,
"T_tank_w_in [°C]": cu.K2C(self.T_tank_w_in_K),
"T_sup_w [°C]": cu.K2C(self.T_sup_w_K),
}
)
if not self.tank_always_full or (self.tank_always_full and self.prevent_simultaneous_flow):
r["tank_level [-]"] = level_solved
return r
# =============================================================
# Template Method Hooks (override in scenario subclasses)
# =============================================================
def _get_activation_flags(
self,
hour_of_day: float,
) -> dict[str, bool]:
"""Return per-subsystem schedule activation flags for *hour_of_day*.
Returns a dict mapping subsystem name → ``True`` if the
subsystem should be active at this hour.
Default: delegates to ``self.stc.is_preheat_on()`` when an
STC is attached (backward-compat); returns ``{}`` otherwise.
Scenario subclasses override to implement custom schedules.
"""
if self.stc is not None:
return {"stc": self.stc.is_preheat_on(hour_of_day)}
return {}
def _needs_solar_input(self) -> bool:
"""Return True if any subsystem requires solar irradiance (I_DN, I_dH).
Default: checks if self.stc or self.pv exists (backward-compat).
Scenario subclasses should override this if they don't attach
components directly to self.stc/self.pv.
"""
return self.stc is not None or self.pv is not None
def _build_residual_fn(
self,
ctx: "StepContext",
ctrl: "ControlState",
dt_s: float,
T_tank_w_in_K_n: float,
T_sup_w_K_n: float,
tank_level: float,
sub_states: dict,
): # -> Callable[[float], float]
"""Return the 1-D energy-balance residual function for *brentq*.
Default implementation: passes *sub_states* as fixed values
(backward-compatible, semi-implicit).
Scenario subclasses override this to re-evaluate their
subsystem physics at ``T_cand`` during every iteration
of the nonlinear solver, achieving a fully implicit solve.
Parameters
----------
ctx : StepContext
Current-step immutable context.
ctrl : ControlState
HP control decisions.
dt_s : float
Time-step size [s].
T_tank_w_in_K_n : float
Mains water inlet temperature [K] (fixed for this step).
T_sup_w_K_n : float
Mains supply temperature [K] (for mixing valve).
tank_level : float
Pre-computed next-step tank level approximation.
sub_states : dict
Subsystem states computed by ``_run_subsystems()``
(frozen at ``T_tank_n``; override to unfreeze).
Returns
-------
Callable[[float], float]
``f(T_cand_K) -> residual`` for use with ``root_scalar``.
"""
def residual(T_cand_K: float) -> float:
return tank_mass_energy_residual(
[T_cand_K, tank_level],
ctx,
ctrl,
dt_s,
T_tank_w_in_K_n,
T_sup_w_K_n,
self.T_mix_w_out_K,
self.C_tank,
self.UA_tank,
self.V_tank_full,
self._subsystems,
sub_states,
)[0]
return residual
def _run_subsystems(
self,
ctx: "StepContext",
ctrl: "ControlState",
dt: float,
T_tank_w_in_K: float,
) -> dict[str, dict]:
"""Step all attached subsystems and return their state dicts.
Default: iterates ``self._subsystems`` (backward-compat).
Scenario subclasses override this to call specific subsystems
without touching ``self._subsystems``.
"""
sub_states: dict[str, dict] = {}
for name, sub in self._subsystems.items():
sub_states[name] = sub.step(ctx, ctrl, dt, T_tank_w_in_K)
return sub_states
def _augment_results(
self,
r: dict,
ctx: "StepContext",
ctrl: "ControlState",
sub_states: dict[str, dict],
T_solved_K: float,
) -> dict:
"""Append subsystem result columns to the step result dict.
Default: iterates ``self._subsystems.assemble_results()`` (backward-compat).
Scenario subclasses override to call specific subsystem assemblers.
"""
for name, sub in self._subsystems.items():
r.update(sub.assemble_results(ctx, ctrl, sub_states[name], T_solved_K))
return r
def _postprocess(self, df: "pd.DataFrame") -> "pd.DataFrame":
"""Post-process the result DataFrame (exergy calculations).
Default: delegates to ``self.postprocess_exergy()`` (backward-compat).
Scenario subclasses override to append subsystem-specific exergy columns
after calling ``super()._postprocess(df)``.
"""
return self.postprocess_exergy(df)
# =============================================================
# Main dynamic simulation
# =============================================================
[docs]
def analyze_dynamic(
self,
simulation_period_sec: int,
dt_s: int,
T_tank_w_init_C: float,
dhw_usage_schedule,
T0_schedule,
I_DN_schedule=None,
I_dH_schedule=None,
T_sup_w_schedule=None,
tank_level_init: float = 1.0,
result_save_csv_path: str | None = None,
) -> pd.DataFrame:
"""Run a time-stepping dynamic simulation.
Fully implicit scheme: ``fsolve`` solves for
``[T_next, level_next]`` each timestep.
Parameters
----------
simulation_period_sec : int
Total simulation duration [s].
dt_s : int
Time step size [s].
T_tank_w_init_C : float
Initial tank temperature [°C].
dhw_usage_schedule : np.ndarray
DHW volumetric flow rate per step [m³/s].
T0_schedule : array-like
Outdoor temperature per step [°C].
I_DN_schedule : array-like | None
Direct-normal irradiance per step [W/m²].
I_dH_schedule : array-like | None
Diffuse-horizontal irradiance [W/m²].
T_sup_w_schedule : array-like | None
Mains water supply temperature per step [°C].
If ``None``, the constructor value ``T_sup_w``
is used for every step (backward-compatible).
tank_level_init : float
Initial fractional tank level (0–1).
result_save_csv_path : str | None
Optional CSV output path.
Returns
-------
pd.DataFrame
Per-timestep result DataFrame.
"""
time: np.ndarray = np.arange(
0,
simulation_period_sec,
dt_s,
)
tN: int = len(time)
T0_schedule = np.array(T0_schedule)
if len(T0_schedule) != tN:
raise ValueError(
f"T0_schedule length ({len(T0_schedule)}) != time length ({tN})",
)
if I_DN_schedule is not None and len(I_DN_schedule) != tN:
raise ValueError(
f"I_DN_schedule length ({len(I_DN_schedule)}) != tN ({tN})",
)
if I_dH_schedule is not None and len(I_dH_schedule) != tN:
raise ValueError(
f"I_dH_schedule length ({len(I_dH_schedule)}) != tN ({tN})",
)
# T_sup_w schedule: fallback to constructor constant
if T_sup_w_schedule is not None:
T_sup_w_arr: np.ndarray = np.array(
T_sup_w_schedule,
dtype=float,
)
if len(T_sup_w_arr) != tN:
raise ValueError(
f"T_sup_w_schedule length ({len(T_sup_w_arr)}) != tN ({tN})",
)
else:
T_sup_w_arr = np.full(tN, self.T_sup_w)
self.time: np.ndarray = time
self.dt: int = dt_s
# DHW schedule handling: direct m³/s flow array
self.dhw_flow_m3s: np.ndarray = np.asarray(
dhw_usage_schedule,
dtype=float,
)
if len(self.dhw_flow_m3s) != tN:
raise ValueError(
f"dhw_usage_schedule length ({len(self.dhw_flow_m3s)}) != tN ({tN})",
)
T_tank_w_K: float = cu.C2K(T_tank_w_init_C)
tank_level: float = tank_level_init
is_refilling: bool = False
hp_is_on_prev: bool = False
results_data: list[dict] = []
# STC/PV schedule flags — kept for StepContext.I_DN/I_dH defaults
_use_solar: bool = self._needs_solar_input()
for n in tqdm(range(tN), desc="ASHPB Simulating"):
t_s: float = time[n]
hr: float = t_s * cu.s2h
hour_of_day: float = (t_s % (24 * cu.h2s)) * cu.s2h
# Per-step mains water supply temperature
T_sup_w_n: float = T_sup_w_arr[n]
T_sup_w_K_n: float = cu.C2K(T_sup_w_n)
T_tank_w_in_K_n: float = T_sup_w_K_n
# Sync self fields for _calc_state compat
self.T_sup_w = T_sup_w_n
self.T_sup_w_K = T_sup_w_K_n
self.T_tank_w_in = T_sup_w_n
self.T_tank_w_in_K = T_tank_w_in_K_n
# Subsystem activation schedule — delegated to Hook
activation_flags: dict[str, bool] = self._get_activation_flags(hour_of_day)
ctx: StepContext = StepContext(
n=n,
current_time_s=t_s,
current_hour=hr,
hour_of_day=hour_of_day,
T0=T0_schedule[n],
T0_K=cu.C2K(T0_schedule[n]),
activation_flags=activation_flags,
T_tank_w_K=T_tank_w_K,
tank_level=tank_level,
dV_mix_w_out=self.dhw_flow_m3s[n],
I_DN=(I_DN_schedule[n] if _use_solar and I_DN_schedule is not None else 0.0),
I_dH=(I_dH_schedule[n] if _use_solar and I_dH_schedule is not None else 0.0),
T_sup_w_K=T_sup_w_K_n,
)
# --- Phase A: control decisions ---
hp_is_on, hp_result, Q_ref_cond = self._determine_hp_state(ctx, hp_is_on_prev)
hp_is_on_prev = hp_is_on
dV_tank_w_in_ctrl, is_refilling = determine_tank_refill_flow(
dt=dt_s,
tank_level=ctx.tank_level,
dV_tank_w_out=self.dV_tank_w_out,
V_tank_full=self.V_tank_full,
tank_always_full=self.tank_always_full,
prevent_simultaneous_flow=self.prevent_simultaneous_flow,
tank_level_lower_bound=self.tank_level_lower_bound,
tank_level_upper_bound=self.tank_level_upper_bound,
dV_tank_w_in_refill=self.dV_tank_w_in_refill,
is_refilling=is_refilling,
)
ctrl: ControlState = ControlState(
is_on=hp_is_on,
Q_heat_source=Q_ref_cond,
dV_tank_w_in_ctrl=dV_tank_w_in_ctrl,
result=hp_result,
)
# --- Phase A-2: subsystem step (via Hook) ---
sub_states: dict[str, dict] = self._run_subsystems(
ctx,
ctrl,
dt_s,
T_tank_w_in_K_n,
)
# --- Phase B: implicit solve (1D over T_next since mass is explicit) ---
# Uncouple mass explicitly:
alp_prev: float = min(
1.0, max(0.0, (self.T_mix_w_out_K - T_sup_w_K_n) / max(1e-6, ctx.T_tank_w_K - T_sup_w_K_n))
)
dV_tank_w_out_prev = alp_prev * ctx.dV_mix_w_out
dV_tank_w_in_prev = dV_tank_w_out_prev if ctrl.dV_tank_w_in_ctrl is None else ctrl.dV_tank_w_in_ctrl
level_next_approx = ctx.tank_level + (dV_tank_w_in_prev - dV_tank_w_out_prev) * dt_s / self.V_tank_full
tank_level = max(0.001, min(1.0, level_next_approx))
residual_1d = self._build_residual_fn(
ctx,
ctrl,
dt_s,
T_tank_w_in_K_n,
T_sup_w_K_n,
tank_level,
sub_states,
)
from scipy.optimize import root_scalar
try:
res = root_scalar(residual_1d, bracket=[cu.C2K(0.0), cu.C2K(100.0)], method="brentq")
converged = getattr(res, "converged", False)
root_val = getattr(res, "root", np.nan)
if converged and not np.isnan(root_val):
T_tank_w_K = float(root_val)
ier = 1
else:
raise ValueError(f"Not converged or NaN: {res}")
except Exception: # Fallback to explicit step if anything fails
# Exception ignored; explicit Euler fallback will correctly handle the state
# Explicit Euler step for energy:
# r_energy = C_curr * T_next - C_curr * T_curr - dt * (Q_total - UA*(T_curr - T0)) = 0
Q_hp_val = ctrl.Q_heat_source
alp_curr = min(
1.0, max(0.0, (self.T_mix_w_out_K - T_sup_w_K_n) / max(1e-6, ctx.T_tank_w_K - T_sup_w_K_n))
)
dV_out_curr = alp_curr * ctx.dV_mix_w_out
Q_flow_curr = c_w * rho_w * dV_out_curr * (T_sup_w_K_n - ctx.T_tank_w_K)
Q_loss_curr = self.UA_tank * (ctx.T_tank_w_K - ctx.T0_K)
Q_tot = Q_hp_val + Q_flow_curr - Q_loss_curr # Assumes sub_total = 0 explicitly for fallback
T_tank_w_K = ctx.T_tank_w_K + dt_s * Q_tot / self.C_tank
ier = 0
if np.isnan(T_tank_w_K) and n < 10:
pass # Silenced NaN fallback debug print
# --- Phase C: core + subsystem results (via Hook) ---
r: dict = self._assemble_core_results(
ctx,
ctrl,
T_tank_w_K,
tank_level,
ier,
)
r = self._augment_results(r, ctx, ctrl, sub_states, T_tank_w_K)
results_data.append(r)
results_df: pd.DataFrame = pd.DataFrame(results_data)
results_df = self._postprocess(results_df)
if result_save_csv_path:
results_df.to_csv(
result_save_csv_path,
index=False,
)
return results_df
# =============================================================
# Exergy post-processing (ASHP-specific)
# =============================================================
[docs]
def postprocess_exergy(self, df: pd.DataFrame) -> pd.DataFrame:
"""Compute ASHP-specific exergy variables.
Owns the full HP exergy topology:
1. Refrigerant state-point exergy (CoolProp)
2. Electricity = exergy (compressor, OU fan, UV)
3. Air exergy (outdoor unit)
4. Heat exchanger Carnot exergy (condenser, evaporator)
5. Water exergy (tank inlet/outlet, mixing valve)
6. Heat loss exergy, tank stored exergy
7. Subsystem exergy via ``calc_exergy()`` protocol
8. Component-level exergy destruction
9. Exergetic efficiency metrics
Parameters
----------
df : pd.DataFrame
Result DataFrame from ``analyze_dynamic()``.
Returns
-------
pd.DataFrame
DataFrame with exergy columns appended.
"""
from .thermodynamics import (
calc_exergy_flow,
calc_refrigerant_exergy,
convert_electricity_to_exergy,
)
df = df.copy()
T0_K = cu.C2K(df["T0 [°C]"])
T_tank_K = cu.C2K(df["T_tank_w [°C]"])
# ── 1. Refrigerant exergy (uses pre-computed h/s from calc_ref_state)
df = calc_refrigerant_exergy(df, self.ref, T0_K)
# ── 2. Electricity = exergy ────────────────────────
df = convert_electricity_to_exergy(df)
# ── 3. Air exergy (outdoor unit) ───────────────────
if "dV_ou_a [m3/s]" in df.columns and "T_ou_a_in [°C]" in df.columns:
G_a = c_a * rho_a * df["dV_ou_a [m3/s]"]
Tin = cu.C2K(df["T_ou_a_in [°C]"])
Tmid = cu.C2K(df["T_ou_a_mid [°C]"])
Tout = cu.C2K(df["T_ou_a_out [°C]"]) if "T_ou_a_out [°C]" in df.columns else Tin
df["X_a_ou_in [W]"] = calc_exergy_flow(G_a, Tin, T0_K)
df["X_a_ou_out [W]"] = calc_exergy_flow(G_a, Tout, T0_K)
df["X_a_ou_mid [W]"] = calc_exergy_flow(G_a, Tmid, T0_K)
# ── 4. HX exergy (Carnot form) ─────────────────────
df["X_ref_cond [W]"] = df["Q_ref_cond [W]"] * (1 - T0_K / cu.C2K(df["T_ref_cond_sat_v [°C]"]))
df["X_ref_evap [W]"] = df["Q_ref_evap [W]"] * (1 - T0_K / cu.C2K(df["T_ref_evap_sat [°C]"]))
# ── 5. Water exergy (inlet / outlet) ───────────────
df["X_tank_w_in [W]"] = calc_exergy_flow(
c_w * rho_w * df["dV_tank_w_in [m3/s]"].fillna(0),
cu.C2K(df["T_tank_w_in [°C]"]),
T0_K,
)
df["X_tank_w_out [W]"] = calc_exergy_flow(
c_w * rho_w * df["dV_tank_w_out [m3/s]"].fillna(0),
cu.C2K(df["T_tank_w [°C]"]),
T0_K,
)
df["X_mix_w_out [W]"] = calc_exergy_flow(
c_w * rho_w * df["dV_mix_w_out [m3/s]"].fillna(0),
cu.C2K(df["T_mix_w_out [°C]"]),
T0_K,
)
df["X_mix_sup_w_in [W]"] = calc_exergy_flow(
c_w * rho_w * df["dV_mix_sup_w_in [m3/s]"].fillna(0),
cu.C2K(df["T_sup_w [°C]"]),
T0_K,
)
# ── 6. Heat loss exergy ────────────────────────────
df["X_tank_loss [W]"] = df["Q_tank_loss [W]"] * (1 - T0_K / T_tank_K)
# ── 7. Tank stored exergy ──────────────────────────
tank_level = df["tank_level [-]"] if "tank_level [-]" in df.columns else 1.0
C_tank_actual = self.C_tank * tank_level
T_tank_K_prev = T_tank_K.shift(1)
df["Xst_tank [W]"] = (1 - T0_K / T_tank_K) * C_tank_actual * (T_tank_K - T_tank_K_prev) / self.dt
df.loc[df.index[0], "Xst_tank [W]"] = 0.0
# ── 8. Removed Subsystem exergy (protocol) ─────────
# Subsystems handle their own exergy via _postprocess hook.
# ── 9. Total exergy input (system-level) ──────────
X_tot = df["E_cmp [W]"] + df["E_ou_fan [W]"]
if "X_uv [W]" in df.columns:
X_tot = X_tot + df["X_uv [W]"].fillna(0)
df["X_tot [W]"] = X_tot
# ── 10. Component exergy destruction ───────────────
# Xc = ΣX_in − ΣX_out ≥ 0 (2nd law)
df["Xc_cmp [W]"] = df["X_cmp [W]"] + df["X_ref_cmp_in [W]"] - df["X_ref_cmp_out [W]"]
df["Xc_ref_cond [W]"] = df["X_ref_cmp_out [W]"] - df["X_ref_exp_in [W]"] - df["X_ref_cond [W]"]
df["Xc_exp [W]"] = df["X_ref_exp_in [W]"] - df["X_ref_exp_out [W]"]
df["Xc_ref_evap [W]"] = (df["X_ref_exp_out [W]"] + df["X_a_ou_in [W]"]) - (
df["X_ref_cmp_in [W]"] + df["X_a_ou_mid [W]"]
)
df["Xc_ou_fan [W]"] = df["X_ou_fan [W]"] + df["X_a_ou_mid [W]"] - df["X_a_ou_out [W]"]
df["Xc_mix [W]"] = df["X_tank_w_out [W]"] + df["X_mix_sup_w_in [W]"] - df["X_mix_w_out [W]"]
# 10g. Storage tank
X_in_tank = df["X_ref_cond [W]"] + df["X_tank_w_in [W]"].fillna(0)
if "X_uv [W]" in df.columns:
X_in_tank = X_in_tank + df["X_uv [W]"].fillna(0)
X_out_tank = df[
"Xst_tank [W]"
] # df['X_tank_loss [W]']를 제외하는 이유는 X_tank_loss 또한 exergy consumption에 포함시기 위함임
if "X_tank_w_out [W]" in df.columns:
X_out_tank = X_out_tank + df["X_tank_w_out [W]"].fillna(0)
df["Xc_tank [W]"] = X_in_tank - X_out_tank
# ── 11. Exergetic efficiency metrics ───────────────
df["X_eff_ref [-]"] = df["X_ref_cond [W]"] / df["X_cmp [W]"].replace(0, np.nan)
df["X_eff_sys [-]"] = df["X_ref_cond [W]"] / df["X_tot [W]"].replace(0, np.nan)
df["X_eff_tank [-]"] = 1 - df["Xc_tank [W]"] / X_in_tank.replace(0, np.nan)
X_in_mix = df["X_tank_w_out [W]"].fillna(0) + df["X_mix_sup_w_in [W]"].fillna(0)
df["X_eff_mix [-]"] = 1 - df["Xc_mix [W]"] / X_in_mix.replace(0, np.nan)
X_in_cmp = df["X_cmp [W]"] + df["X_ref_cmp_in [W]"]
df["X_eff_cmp [-]"] = 1 - df["Xc_cmp [W]"] / X_in_cmp.replace(0, np.nan)
df["X_eff_ref_cond [-]"] = 1 - df["Xc_ref_cond [W]"] / df["X_ref_cmp_out [W]"].replace(0, np.nan)
df["X_eff_exp [-]"] = 1 - df["Xc_exp [W]"] / df["X_ref_exp_in [W]"].replace(0, np.nan)
a_ou_in = df["X_a_ou_in [W]"].fillna(0) if "X_a_ou_in [W]" in df.columns else 0.0
X_in_ref_evap = df["X_ref_exp_out [W]"] + a_ou_in
df["X_eff_ref_evap [-]"] = 1 - df["Xc_ref_evap [W]"] / X_in_ref_evap.replace(0, np.nan)
a_ou_mid = df["X_a_ou_mid [W]"].fillna(0) if "X_a_ou_mid [W]" in df.columns else 0.0
X_in_ou_fan = df["X_ou_fan [W]"] + a_ou_mid
df["X_eff_ou_fan [-]"] = 1 - df["Xc_ou_fan [W]"] / X_in_ou_fan.replace(0, np.nan)
return df