diff --git a/examples/35_etes/run_etes_performance.py b/examples/35_etes/run_etes_performance.py new file mode 100644 index 000000000..fb1e55ee1 --- /dev/null +++ b/examples/35_etes/run_etes_performance.py @@ -0,0 +1,198 @@ +""" +Example showcasing the ETESPerformanceModel. + +This example simulates one week (168 hours) of an Electric Thermal Energy +Storage (ETES) system serving a steady industrial thermal load. The +charging electricity profile mimics a renewables-heavy grid: ample, cheap +electricity at night and midday (solar/wind), with a deficit during the +evening peak. The script runs both P-ETES (decoupled) and R-ETES +(integrated) configurations and produces a plot of state-of-charge, +charging / discharging rates, and heat delivered vs. demand. + +Run with the h2integrate conda environment: + + conda run -n h2integrate python run_etes_performance.py +""" + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import openmdao.api as om + +from h2integrate.converters.heat.etes import ETESPerformanceModel + + +# --------------------------------------------------------------------------- +# Scenario setup +# --------------------------------------------------------------------------- +N_HOURS = 168 # 1 week +DT_S = 3600 # 1-hour timestep +HEAT_LOAD_KW = 10_000.0 # steady 10 MW_th industrial load + + +def build_electricity_profile(n: int) -> np.ndarray: + """Synthetic available-electricity profile (kW_e) with diurnal cycles. + + - High at night (cheap wind) + - Peak around midday (solar) + - Low during evening ramp (no sun, low wind) + """ + t = np.arange(n) + hour_of_day = t % 24 + # Solar bump (Gaussian around hour 13) + solar = 60_000.0 * np.exp(-0.5 * ((hour_of_day - 13) / 3.5) ** 2) + # Wind: noisy baseline with a dip in the evening (hours 17-21) + rng = np.random.default_rng(42) + wind = 25_000.0 + 5_000.0 * np.sin(t / 7.0) + rng.normal(0, 2_000.0, n) + evening_dip = np.where((hour_of_day >= 17) & (hour_of_day <= 21), -15_000.0, 0.0) + profile = np.clip(solar + wind + evening_dip, 0.0, None) + return profile + + +# --------------------------------------------------------------------------- +# Model configurations (parameter values from MILP ETES doc) +# --------------------------------------------------------------------------- +PETES_PARAMS = { + "etes_type": "P-ETES", + "S_TES_kWh": 600_000.0, # 600 MWh thermal storage + "S_ch_kW": 100_000.0, # 100 MW electric charging unit + "S_dis_kW": 20_000.0, # 20 MW thermal discharging unit + "eta_ch": 0.95, # ENDURING particle heater + "eta_dis": 0.73, # MPBHX + "f_loss": 0.0068, # per-hour fractional loss + "SOC_min": 0.05, + "SOC_init": 0.5, + "cost_year": 2026, +} + +RETES_PARAMS = { + "etes_type": "R-ETES", + "S_TES_kWh": 600_000.0, + "eta_ch": 0.98, # refractory heater + "eta_dis": 0.90, + "f_loss": 0.0068, + "f_ch_max": 0.3, + "f_dis_max": 0.2, + "SOC_min": 0.05, + "SOC_init": 0.5, + "cost_year": 2026, +} + + +def run_etes(perf_params: dict, electricity_in_kW: np.ndarray, heat_demand_kW: np.ndarray): + """Build an OpenMDAO problem, run the ETES performance model, return the + populated problem instance.""" + plant_config = { + "plant": { + "simulation": {"n_timesteps": N_HOURS, "dt": DT_S}, + } + } + tech_config = {"model_inputs": {"performance_parameters": perf_params}} + + prob = om.Problem() + prob.model.add_subsystem( + "etes", + ETESPerformanceModel( + tech_config=tech_config, + plant_config=plant_config, + driver_config={}, + ), + promotes=["*"], + ) + prob.setup() + prob.set_val("electricity_in_kW", electricity_in_kW) + prob.set_val("heat_demand_kW", heat_demand_kW) + prob.run_model() + return prob + + +def summarize(label: str, prob: om.Problem): + heat_out = prob.get_val("heat_out_kW", units="kW") + unmet = prob.get_val("unmet_heat_demand_kW", units="kW") + elec_used = prob.get_val("electricity_consumed_kW", units="kW") + rte = float(prob.get_val("round_trip_efficiency")[0]) + total_delivered = float(prob.get_val("E_total_delivered_kWh", units="kW*h")[0]) + + print(f"\n--- {label} ---") + print(f" Total heat delivered : {total_delivered/1e3:,.1f} MWh_th") + print(f" Total unmet demand : {np.sum(unmet)/1e3:,.1f} MWh_th " + f"({np.sum(unmet) / np.sum(heat_out + unmet) * 100:.2f}% of demand)") + print(f" Total electricity in : {np.sum(elec_used)/1e3:,.1f} MWh_e") + print(f" Round-trip efficiency: {rte * 100:.1f}%") + + +def plot_results(electricity_in: np.ndarray, heat_demand: np.ndarray, + petes_prob: om.Problem, retes_prob: om.Problem, + outfile: Path): + fig, axes = plt.subplots(4, 1, figsize=(12, 11), sharex=True) + t = np.arange(N_HOURS) + + # 1. Available electricity + axes[0].plot(t, electricity_in / 1e3, color="tab:blue") + axes[0].set_ylabel("Available\nelectricity (MW_e)") + axes[0].set_title("Scenario: 1-week renewables-heavy electricity supply, 10 MW_th steady load") + axes[0].grid(True, alpha=0.3) + + # 2. SOC + for label, prob, color in [ + ("P-ETES", petes_prob, "tab:orange"), + ("R-ETES", retes_prob, "tab:green"), + ]: + E_st = prob.get_val("E_st_kWh", units="kW*h") + S_TES = (PETES_PARAMS if label == "P-ETES" else RETES_PARAMS)["S_TES_kWh"] + axes[1].plot(t, E_st / S_TES * 100, label=label, color=color) + axes[1].set_ylabel("State of charge (%)") + axes[1].legend(loc="upper right") + axes[1].grid(True, alpha=0.3) + axes[1].set_ylim(0, 105) + + # 3. Charging / discharging rates + for label, prob, color in [ + ("P-ETES", petes_prob, "tab:orange"), + ("R-ETES", retes_prob, "tab:green"), + ]: + Q_ch = prob.get_val("Q_ch_kW", units="kW") + Q_dis = prob.get_val("Q_dis_kW", units="kW") + axes[2].plot(t, Q_ch / 1e3, label=f"{label} charge", color=color, linestyle="-") + axes[2].plot(t, -Q_dis / 1e3, label=f"{label} discharge", color=color, linestyle="--") + axes[2].axhline(0, color="black", linewidth=0.5) + axes[2].set_ylabel("Storage flow\n(MW_th, +charge / -discharge)") + axes[2].legend(loc="upper right", ncol=2, fontsize=8) + axes[2].grid(True, alpha=0.3) + + # 4. Heat delivered vs. demand + axes[3].plot(t, heat_demand / 1e3, color="black", linewidth=2, label="Demand") + for label, prob, color in [ + ("P-ETES", petes_prob, "tab:orange"), + ("R-ETES", retes_prob, "tab:green"), + ]: + heat_out = prob.get_val("heat_out_kW", units="kW") + axes[3].plot(t, heat_out / 1e3, label=f"{label} delivered", + color=color, linestyle="--", alpha=0.85) + axes[3].set_ylabel("Thermal power (MW_th)") + axes[3].set_xlabel("Time (h)") + axes[3].legend(loc="lower right", fontsize=8) + axes[3].grid(True, alpha=0.3) + + fig.tight_layout() + fig.savefig(outfile, dpi=120) + print(f"\nWrote plot: {outfile}") + + +def main(): + electricity_in = build_electricity_profile(N_HOURS) + heat_demand = np.full(N_HOURS, HEAT_LOAD_KW) + + petes_prob = run_etes(PETES_PARAMS, electricity_in, heat_demand) + retes_prob = run_etes(RETES_PARAMS, electricity_in, heat_demand) + + summarize("P-ETES (decoupled)", petes_prob) + summarize("R-ETES (integrated)", retes_prob) + + out = Path(__file__).parent / "etes_performance_results.png" + plot_results(electricity_in, heat_demand, petes_prob, retes_prob, out) + + +if __name__ == "__main__": + main() diff --git a/h2integrate/converters/heat/__init__.py b/h2integrate/converters/heat/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/h2integrate/converters/heat/etes.py b/h2integrate/converters/heat/etes.py new file mode 100644 index 000000000..3d715c2d4 --- /dev/null +++ b/h2integrate/converters/heat/etes.py @@ -0,0 +1,271 @@ +""" +Electric Thermal Energy Storage (ETES) performance model. + +Implements the energy-balance equations from the MILP ETES Optimization +Formulation (Lidor, 2026). Two ETES types are supported: + +- "P-ETES": Particle-type ETES with decoupled charging unit (electric heater), + storage capacity, and discharging unit (e.g., MPBHX). Charging and + discharging rates are limited by their respective unit power ratings + (S_ch, S_dis). +- "R-ETES": Refractory-type ETES (e.g., Rondo, Antora) where charging and + discharging rates are coupled to the storage capacity via rate constants + f_ch_max and f_dis_max. + +The dispatch is implemented as a simple heuristic (not the MILP): +at each timestep, the model attempts to meet the thermal load by discharging +from storage, then uses any available electricity to charge the storage, all +subject to the capacity, rate, and SOC constraints from the formulation. +""" + +import numpy as np +import openmdao.api as om +from attrs import field, define + +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.core.validators import gt_zero, gte_zero + + +def _validate_etes_type(instance, attribute, value): + valid = ("P-ETES", "R-ETES") + if value not in valid: + raise ValueError(f"etes_type must be one of {valid}, got {value!r}") + + +@define(kw_only=True) +class ETESPerformanceModelConfig(BaseConfig): + """ + Configuration class for the ETESPerformanceModel. + + Args: + etes_type (str): ETES type. One of "P-ETES" (decoupled charging, + storage, and discharging components) or "R-ETES" (integrated + with rate-coupled charging/discharging). + S_TES_kWh (float): ETES thermal storage capacity in kWh_th. + S_ch_kW (float): Charging unit electric power rating in kW_e + (used for P-ETES only). + S_dis_kW (float): Discharging unit thermal power rating in kW_th + (used for P-ETES only). + eta_ch (float): Charging efficiency (0 < eta_ch <= 1). + eta_dis (float): Discharging efficiency (0 < eta_dis <= 1). + f_loss (float): Per-timestep storage loss fraction (fraction of + stored energy lost each timestep). + SOC_min (float): Minimum state of charge (fraction of S_TES). + SOC_init (float): Initial state of charge (fraction of S_TES) at + the beginning of the simulation. Defaults to 0.5 per the doc. + f_ch_max (float): Maximum charging rate as a fraction of S_TES per + hour (used for R-ETES only). + f_dis_max (float): Maximum discharging rate as a fraction of S_TES + per hour (used for R-ETES only). + cost_year (int): Year for cost estimation. + """ + + etes_type: str = field(validator=_validate_etes_type) + S_TES_kWh: float = field(validator=gt_zero) + eta_ch: float = field(validator=gt_zero) + eta_dis: float = field(validator=gt_zero) + f_loss: float = field(validator=gte_zero) + SOC_min: float = field(default=0.0, validator=gte_zero) + SOC_init: float = field(default=0.5, validator=gte_zero) + S_ch_kW: float = field(default=0.0, validator=gte_zero) + S_dis_kW: float = field(default=0.0, validator=gte_zero) + f_ch_max: float = field(default=0.0, validator=gte_zero) + f_dis_max: float = field(default=0.0, validator=gte_zero) + cost_year: int = field(validator=gt_zero) + + +class ETESPerformanceModel(om.ExplicitComponent): + """ + An OpenMDAO component that simulates ETES dispatch over a time-series + horizon, given an available electricity profile (for charging) and a + thermal load profile. + + Inputs (time-series): + electricity_in_kW: Electricity available for charging the ETES at + each timestep, in kW_e. + heat_demand_kW: Thermal load demand at each timestep, in kW_th. + + Outputs (time-series): + heat_out_kW: Thermal energy delivered to load at each timestep, in + kW_th. + unmet_heat_demand_kW: Unmet thermal load at each timestep, in + kW_th. + electricity_consumed_kW: Electricity actually consumed for + charging at each timestep, in kW_e. + E_st_kWh: Stored thermal energy at each timestep, in kWh_th. + Q_ch_kW, Q_dis_kW: Charging / discharging rates (thermal), in + kW_th. + Q_st_loss_kW, Q_ch_loss_kW, Q_dis_loss_kW: Storage / charging / + discharging loss rates (thermal), in kW_th. + + Scalar outputs: + E_total_delivered_kWh: Total thermal energy delivered to load + over the simulation, in kWh_th. + overall_round_trip_efficiency: Ratio of total thermal energy delivered to + total electricity consumed. + """ + + # TODO: Change all thermal energies to Q and reserve E for electricity + + def initialize(self): + self.options.declare("driver_config", types=dict) + self.options.declare("plant_config", types=dict) + self.options.declare("tech_config", types=dict) + + def setup(self): + self.config = ETESPerformanceModelConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), + additional_cls_name=self.__class__.__name__, + ) + + sim = self.options["plant_config"]["plant"]["simulation"] + self.n_timesteps = int(sim["n_timesteps"]) + self.dt_h = float(sim["dt"]) / 3600.0 # timestep in hours + + # Time-series inputs + self.add_input( + "electricity_in_kW", + val=0.0, + shape=self.n_timesteps, + units="kW", + desc="Electricity available for charging at each timestep", + ) + self.add_input( + "heat_demand_kW", + val=0.0, + shape=self.n_timesteps, + units="kW", + desc="Thermal load demand at each timestep", + ) + + # Time-series outputs + for name, units, desc in [ + ("heat_out_kW", "kW", "Thermal energy delivered to load"), + ("unmet_heat_demand_kW", "kW", "Unmet thermal load"), + ("electricity_consumed_kW", "kW", "Electricity consumed for charging"), + ("E_st_kWh", "kW*h", "Stored thermal energy at end of each timestep"), + ("Q_ch_kW", "kW", "Thermal charging rate"), + ("Q_dis_kW", "kW", "Thermal discharging rate (withdrawn from storage)"), + ("Q_st_loss_kW", "kW", "Storage loss rate"), + ("Q_ch_loss_kW", "kW", "Charging loss rate"), + ("Q_dis_loss_kW", "kW", "Discharging loss rate"), + ]: + self.add_output(name, val=0.0, shape=self.n_timesteps, units=units, desc=desc) + + # Scalar outputs + self.add_output( + "E_total_delivered_kWh", + val=0.0, + units="kW*h", + desc="Total thermal energy delivered to load over the simulation", + ) + self.add_output( + "overall_round_trip_efficiency", + val=0.0, + desc="Ratio of thermal energy delivered to electricity consumed", + ) + + def compute(self, inputs, outputs): + cfg = self.config + n = self.n_timesteps + dt = self.dt_h + + S_TES = float(cfg.S_TES_kWh) + eta_ch = float(cfg.eta_ch) + eta_dis = float(cfg.eta_dis) + f_loss = float(cfg.f_loss) + SOC_min = float(cfg.SOC_min) + + E_min = SOC_min * S_TES + E_max = S_TES + + # Rate caps (kW_th withdrawn from / added to storage) + if cfg.etes_type == "P-ETES": + # P-ETES: decoupled. Charging rate cap = S_ch * eta_ch (thermal added to + # storage from electric input of S_ch). Discharging rate cap = S_dis / + # eta_dis (storage withdrawal to deliver S_dis of thermal output). + ch_rate_cap = float(cfg.S_ch_kW) * eta_ch + dis_rate_cap = float(cfg.S_dis_kW) / eta_dis if eta_dis > 0 else np.inf + else: # R-ETES: integrated, rates coupled to S_TES + ch_rate_cap = float(cfg.f_ch_max) * S_TES + dis_rate_cap = float(cfg.f_dis_max) * S_TES + + elec_in = np.asarray(inputs["electricity_in_kW"], dtype=float) + load = np.asarray(inputs["heat_demand_kW"], dtype=float) + + heat_out = np.zeros(n) + unmet = np.zeros(n) + elec_used = np.zeros(n) + E_st = np.zeros(n) + Q_ch = np.zeros(n) + Q_dis = np.zeros(n) + Q_st_loss = np.zeros(n) + Q_ch_loss = np.zeros(n) + Q_dis_loss = np.zeros(n) + + E_prev = float(cfg.SOC_init) * S_TES + + for t in range(n): + # Storage standing loss (energy lost during this timestep, based on + # state at start of timestep). Q_st_loss is the rate over Δt. + loss_energy = E_prev * f_loss + q_st_loss_rate = loss_energy / dt if dt > 0 else 0.0 + + # --- Discharge to meet thermal load --- + # Required storage-side withdrawal rate to meet load: + dis_req_rate = load[t] / eta_dis if eta_dis > 0 else 0.0 # kW_th + # Available energy in storage after standing loss, above SOC_min: + available_E = max(E_prev - loss_energy - E_min, 0.0) + dis_available_rate = available_E / dt if dt > 0 else 0.0 + dis_rate = min(dis_req_rate, dis_rate_cap, dis_available_rate) + heat_delivered = dis_rate * eta_dis + + # --- Charge with available electricity --- + # Thermal added to storage per electricity input: + ch_req_rate_thermal = elec_in[t] * eta_ch + # Headroom in storage after standing loss and discharge: + headroom_E = max( + E_max - (E_prev - loss_energy - dis_rate * dt), + 0.0, + ) + ch_headroom_rate = headroom_E / dt if dt > 0 else 0.0 + ch_rate = min(ch_req_rate_thermal, ch_rate_cap, ch_headroom_rate) + elec_consumed = ch_rate / eta_ch if eta_ch > 0 else 0.0 + + # Update storage state + E_new = E_prev - loss_energy + (ch_rate - dis_rate) * dt + # Numerical clamp - TODO: Is this actually necessary? Test with an extreme case + E_new = min(max(E_new, E_min), E_max) + + # Record outputs + heat_out[t] = heat_delivered + unmet[t] = max(load[t] - heat_delivered, 0.0) + elec_used[t] = elec_consumed + E_st[t] = E_new + Q_ch[t] = ch_rate + Q_dis[t] = dis_rate + Q_st_loss[t] = q_st_loss_rate + # Loss rates from formulation: + # Q_ch_loss = Q_ch * (1 - eta_ch) / eta_ch (electric-side loss) + # Q_dis_loss = Q_dis * (1 - eta_dis) (thermal-side loss) + Q_ch_loss[t] = ch_rate * (1.0 - eta_ch) / eta_ch if eta_ch > 0 else 0.0 + Q_dis_loss[t] = dis_rate * (1.0 - eta_dis) + + E_prev = E_new + + outputs["heat_out_kW"] = heat_out + outputs["unmet_heat_demand_kW"] = unmet + outputs["electricity_consumed_kW"] = elec_used + outputs["E_st_kWh"] = E_st + outputs["Q_ch_kW"] = Q_ch + outputs["Q_dis_kW"] = Q_dis + outputs["Q_st_loss_kW"] = Q_st_loss + outputs["Q_ch_loss_kW"] = Q_ch_loss + outputs["Q_dis_loss_kW"] = Q_dis_loss + + total_heat = float(np.sum(heat_out) * dt) + total_elec = float(np.sum(elec_used) * dt) + outputs["E_total_delivered_kWh"] = total_heat + outputs["overall_round_trip_efficiency"] = ( + total_heat / total_elec if total_elec > 0 else 0.0 + ) diff --git a/h2integrate/converters/heat/etes_cost_model.py b/h2integrate/converters/heat/etes_cost_model.py new file mode 100644 index 000000000..75c03a22a --- /dev/null +++ b/h2integrate/converters/heat/etes_cost_model.py @@ -0,0 +1,156 @@ +""" +Cost model for the Electric Thermal Energy Storage (ETES). + +Implements the linearized cost functions from the MILP ETES Optimization +Formulation (Lidor, 2026): + + C_TES_total = max(C_lin_TES * S_TES + C_const_TES, C_min_TES * S_TES) + C_ch_total = C_lin_ch * S_ch + C_const_ch (P-ETES only) + C_dis_total = C_lin_dis * S_dis + C_const_dis (P-ETES only) + +For R-ETES systems the charging/discharging components are integrated with +storage and S_ch / S_dis cost terms are not used (set their linear/constant +coefficients to zero or leave at default). + +OpEx is taken as a fixed fraction of CapEx per the same convention as +``shell_tube_hx_cost_model.py``. +""" + +from attrs import field, define + +from h2integrate.core.utilities import merge_shared_inputs +from h2integrate.core.validators import gte_zero +from h2integrate.core.model_baseclasses import CostModelBaseClass, CostModelBaseConfig + + +@define(kw_only=True) +class ETESCostModelConfig(CostModelBaseConfig): + """ + Configuration class for the ETESCostModel. + + Args: + C_lin_TES (float): ETES storage linear cost coefficient ($/kWh_th). + C_const_TES (float): ETES storage constant cost term ($). + C_min_TES (float): Minimum allowed storage cost per kWh_th + ($/kWh_th). Acts as a lower bound on (C_lin_TES * S_TES + + C_const_TES) / S_TES. + C_lin_ch (float): Charging unit linear cost coefficient ($/kW_e). + Used for P-ETES only. + C_const_ch (float): Charging unit constant cost term ($). + C_lin_dis (float): Discharging unit linear cost coefficient + ($/kW_th). Used for P-ETES only. + C_const_dis (float): Discharging unit constant cost term ($). + opex_fraction (float): Annual OpEx as a fraction of total CapEx. + """ + + C_lin_TES: float = field(validator=gte_zero) + C_const_TES: float = field(default=0.0, validator=gte_zero) + C_min_TES: float = field(default=0.0, validator=gte_zero) + C_lin_ch: float = field(default=0.0, validator=gte_zero) + C_const_ch: float = field(default=0.0, validator=gte_zero) + C_lin_dis: float = field(default=0.0, validator=gte_zero) + C_const_dis: float = field(default=0.0, validator=gte_zero) + opex_fraction: float = field(default=0.04, validator=gte_zero) + + +class ETESCostModel(CostModelBaseClass): + """ + Linearized ETES cost model. + + Inputs: + S_TES_kWh: Storage capacity (kWh_th). + S_ch_kW: Charging unit electric power rating (kW_e). Used for + P-ETES; set to 0 for R-ETES. + S_dis_kW: Discharging unit thermal power rating (kW_th). Used for + P-ETES; set to 0 for R-ETES. + + Outputs: + CapEx: Total capital expenditure (USD). + OpEx: Annual operating expenditure (USD/year). + """ + + def setup(self): + self.config = ETESCostModelConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "cost"), + additional_cls_name=self.__class__.__name__, + ) + super().setup() + + self.add_input("S_TES_kWh", val=0.0, units="kW*h", desc="Storage capacity") + self.add_input("S_ch_kW", val=0.0, units="kW", desc="Charging unit electric power rating") + self.add_input( + "S_dis_kW", val=0.0, units="kW", desc="Discharging unit thermal power rating" + ) + + self.add_input( + "C_lin_TES", + val=self.config.C_lin_TES, + units="USD/(kW*h)", + desc="Storage linear cost coefficient", + ) + self.add_input( + "C_const_TES", + val=self.config.C_const_TES, + units="USD", + desc="Storage constant cost term", + ) + self.add_input( + "C_min_TES", + val=self.config.C_min_TES, + units="USD/(kW*h)", + desc="Minimum storage cost per kWh_th", + ) + self.add_input( + "C_lin_ch", + val=self.config.C_lin_ch, + units="USD/kW", + desc="Charging unit linear cost coefficient", + ) + self.add_input( + "C_const_ch", + val=self.config.C_const_ch, + units="USD", + desc="Charging unit constant cost term", + ) + self.add_input( + "C_lin_dis", + val=self.config.C_lin_dis, + units="USD/kW", + desc="Discharging unit linear cost coefficient", + ) + self.add_input( + "C_const_dis", + val=self.config.C_const_dis, + units="USD", + desc="Discharging unit constant cost term", + ) + + def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): + S_TES = float(inputs["S_TES_kWh"][0]) + S_ch = float(inputs["S_ch_kW"][0]) + S_dis = float(inputs["S_dis_kW"][0]) + + C_lin_TES = float(inputs["C_lin_TES"][0]) + C_const_TES = float(inputs["C_const_TES"][0]) + C_min_TES = float(inputs["C_min_TES"][0]) + C_lin_ch = float(inputs["C_lin_ch"][0]) + C_const_ch = float(inputs["C_const_ch"][0]) + C_lin_dis = float(inputs["C_lin_dis"][0]) + C_const_dis = float(inputs["C_const_dis"][0]) + + # Linearized storage cost with floor (from the MILP formulation: + # C_TES_total - C_min_TES * S_TES + C_const_TES >= 0 => + # C_TES_total >= C_min_TES * S_TES - C_const_TES) + c_tes_linear = C_lin_TES * S_TES + C_const_TES + c_tes_floor = C_min_TES * S_TES + c_tes = max(c_tes_linear, c_tes_floor) + + c_ch = C_lin_ch * S_ch + (C_const_ch if S_ch > 0 else 0.0) + c_dis = C_lin_dis * S_dis + (C_const_dis if S_dis > 0 else 0.0) + + capex = c_tes + c_ch + c_dis + opex = self.config.opex_fraction * capex + + outputs["CapEx"] = capex + outputs["OpEx"] = opex + # TODO: Expand to include cost of electricity used from heuristic etes model diff --git a/h2integrate/converters/heat/etes_milp.py b/h2integrate/converters/heat/etes_milp.py new file mode 100644 index 000000000..83b8bdcd3 --- /dev/null +++ b/h2integrate/converters/heat/etes_milp.py @@ -0,0 +1,467 @@ +""" +MILP optimizer for the Electric Thermal Energy Storage (ETES). + +Implements the optimization formulation from the MILP ETES Optimization +Formulation (Lidor, 2026). Given a time series of electricity prices and +thermal loads, this module solves for the optimal ETES sizing and dispatch +that minimizes the total annualized cost (annualized CapEx + OpEx + grid +electricity cost). + +Two sizing modes are supported via ``etes_type``: + +- "P-ETES": separate charging unit (S_ch), storage capacity (S_TES), and + discharging unit (S_dis) are all sized independently. +- "R-ETES": only the storage capacity (S_TES) is sized; the charging and + discharging rates are coupled to S_TES via rate constants f_ch_max and + f_dis_max. + +The model is technically a linear program as written in the doc (no +binary variables are required). Mutual-exclusion of charge/discharge can +optionally be enforced with binaries, making the problem a true MILP. + +Example +------- +>>> from h2integrate.converters.heat.etes_milp import ETESMILPConfig, solve_etes_milp +>>> import numpy as np +>>> n = 24 +>>> cfg = ETESMILPConfig( +... etes_type="P-ETES", +... eta_ch=0.95, +... eta_dis=0.73, +... f_loss=0.0068, +... t_ch_min_h=4.0, +... C_lin_TES=5.0, +... C_min_TES=1.5, +... C_lin_ch=100.0, +... C_lin_dis=150.0, +... fixed_charge_rate=0.10, +... ) +>>> price = np.array([0.02] * 8 + [0.10] * 8 + [0.05] * 8) # $/kWh_e +>>> load = np.full(24, 10_000.0) # kW_th +>>> result = solve_etes_milp(cfg, price, load, dt_h=1.0) +>>> result.S_TES_kWh # doctest: +SKIP +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +import pyomo.environ as pyo +from attrs import field as attrs_field, define + +from h2integrate.core.utilities import BaseConfig +from h2integrate.core.validators import gt_zero, gte_zero + + +_PREFERRED_SOLVERS = ("appsi_highs", "highs", "cbc", "glpk") + + +def _validate_etes_type(instance, attribute, value): + valid = ("P-ETES", "R-ETES") + if value not in valid: + raise ValueError(f"etes_type must be one of {valid}, got {value!r}") + + +@define(kw_only=True) +class ETESMILPConfig(BaseConfig): + """ + Configuration for the ETES MILP optimizer. + + Physical parameters: + etes_type (str): "P-ETES" or "R-ETES". + eta_ch (float): Charging efficiency. + eta_dis (float): Discharging efficiency. + f_loss (float): Per-timestep storage loss fraction. + SOC_min (float): Minimum state of charge (fraction of S_TES). + SOC_init (float): Initial state of charge (fraction of S_TES). + Defaults to 0.5 per the doc. + # TODO: Don't think we want the cyclic behavior below + cyclic (bool): If True, force E_st at end of horizon equal to the + initial state (typical for annual operation). Defaults to True. + t_ch_min_h (float): Minimum time (h) for fully charging the + P-ETES; constrains S_ch <= S_TES / (eta_ch * t_ch_min_h). + P-ETES only. + f_ch_max (float): Max charging rate as fraction of S_TES (1/h). + R-ETES only. + f_dis_max (float): Max discharging rate as fraction of S_TES + (1/h). R-ETES only. + + Sizing bounds (optional): + S_TES_min_kWh, S_TES_max_kWh: bounds on storage capacity. + S_ch_min_kW, S_ch_max_kW: bounds on charging unit (P-ETES). + S_dis_min_kW, S_dis_max_kW: bounds on discharging unit (P-ETES). + S_TES_fixed_kWh, S_ch_fixed_kW, S_dis_fixed_kW: if set (non-None), + the corresponding size is fixed (dispatch-only mode). + + Cost parameters (linearized): + C_lin_TES (float): Storage linear cost ($/kWh_th). + C_const_TES (float): Storage constant cost term ($). + C_min_TES (float): Minimum allowed cost per kWh_th ($/kWh_th). + C_lin_ch (float): Charging unit linear cost ($/kW_e). P-ETES. + C_const_ch (float): Charging unit constant cost ($). P-ETES. + C_lin_dis (float): Discharging unit linear cost ($/kW_th). P-ETES. + C_const_dis (float): Discharging unit constant cost ($). P-ETES. + + Annualization: + fixed_charge_rate (float): Annual fixed charges per dollar of + CapEx. Defaults to 0.10. + opex_fraction (float): Annual OpEx as fraction of CapEx. Defaults + to 0.04. + + Solver options: + allow_unmet_load (bool): If True, allow load to be partially met + with a penalty cost. Defaults to False. + unmet_load_penalty (float): Penalty cost per kWh_th of unmet load + ($/kWh_th). Only used if allow_unmet_load is True. + enforce_mutex_charge_dis (bool): If True, add binary variables so + charging and discharging cannot occur simultaneously (true + MILP). Defaults to False. + solver (str | None): Pyomo solver name. If None, the first + available solver from {appsi_highs, highs, cbc, glpk} is used. + solver_options (dict | None): Options passed to the solver. + """ + + etes_type: str = attrs_field(validator=_validate_etes_type) + eta_ch: float = attrs_field(validator=gt_zero) + eta_dis: float = attrs_field(validator=gt_zero) + f_loss: float = attrs_field(validator=gte_zero) + SOC_min: float = attrs_field(default=0.0, validator=gte_zero) + SOC_init: float = attrs_field(default=0.5, validator=gte_zero) + cyclic: bool = attrs_field(default=True) + + t_ch_min_h: float = attrs_field(default=0.0, validator=gte_zero) + f_ch_max: float = attrs_field(default=0.0, validator=gte_zero) + f_dis_max: float = attrs_field(default=0.0, validator=gte_zero) + + S_TES_min_kWh: float = attrs_field(default=0.0, validator=gte_zero) + S_TES_max_kWh: float = attrs_field(default=1.0e12, validator=gt_zero) + S_ch_min_kW: float = attrs_field(default=0.0, validator=gte_zero) + S_ch_max_kW: float = attrs_field(default=1.0e9, validator=gt_zero) + S_dis_min_kW: float = attrs_field(default=0.0, validator=gte_zero) + S_dis_max_kW: float = attrs_field(default=1.0e9, validator=gt_zero) + S_TES_fixed_kWh: float | None = attrs_field(default=None) + S_ch_fixed_kW: float | None = attrs_field(default=None) + S_dis_fixed_kW: float | None = attrs_field(default=None) + + C_lin_TES: float = attrs_field(default=0.0, validator=gte_zero) + C_const_TES: float = attrs_field(default=0.0, validator=gte_zero) + C_min_TES: float = attrs_field(default=0.0, validator=gte_zero) + C_lin_ch: float = attrs_field(default=0.0, validator=gte_zero) + C_const_ch: float = attrs_field(default=0.0, validator=gte_zero) + C_lin_dis: float = attrs_field(default=0.0, validator=gte_zero) + C_const_dis: float = attrs_field(default=0.0, validator=gte_zero) + + fixed_charge_rate: float = attrs_field(default=0.10, validator=gte_zero) + opex_fraction: float = attrs_field(default=0.04, validator=gte_zero) + + allow_unmet_load: bool = attrs_field(default=False) + unmet_load_penalty: float = attrs_field(default=1.0e4) + enforce_mutex_charge_dis: bool = attrs_field(default=False) + solver: str | None = attrs_field(default=None) + solver_options: dict | None = attrs_field(default=None) + + +@dataclass +class ETESMILPResult: + """Result returned by :func:`solve_etes_milp`.""" + + # Optimal sizes + S_TES_kWh: float + S_ch_kW: float + S_dis_kW: float + # Time series (length n_timesteps) + E_st_kWh: np.ndarray + Q_ch_kW: np.ndarray + Q_dis_kW: np.ndarray + Q_st_loss_kW: np.ndarray + Q_ch_loss_kW: np.ndarray + Q_dis_loss_kW: np.ndarray + P_grid_kW: np.ndarray + unmet_load_kW: np.ndarray + # Cost breakdown + capex_USD: float + annual_opex_USD: float + annual_electricity_cost_USD: float + annualized_capex_USD: float + total_annualized_cost_USD: float + # Diagnostics + solver_status: str + termination_condition: str + objective_value: float + + +def _pick_solver(name: str | None): + if name is not None: + solver = pyo.SolverFactory(name) + if not solver.available(exception_flag=False): + raise RuntimeError(f"Requested solver {name!r} is not available.") + return solver, name + for s in _PREFERRED_SOLVERS: + try: + solver = pyo.SolverFactory(s) + if solver.available(exception_flag=False): + return solver, s + except (ValueError, AttributeError): + continue + raise RuntimeError( + f"No LP/MILP solver available; tried {_PREFERRED_SOLVERS}. " + "Install one (e.g. `pip install highspy` for HiGHS)." + ) + + +def solve_etes_milp( + config: ETESMILPConfig, + electricity_price: np.ndarray, + heat_load_kW: np.ndarray, + dt_h: float = 1.0, +) -> ETESMILPResult: + """ + Solve the MILP ETES sizing-and-dispatch problem. + + Args: + config: :class:`ETESMILPConfig` with physical, cost, and solver + parameters. + electricity_price: Array of length ``n`` with grid electricity + prices in $/kWh_e for each timestep. + heat_load_kW: Array of length ``n`` with thermal load demand in + kW_th for each timestep. + dt_h: Timestep length in hours. Defaults to 1.0. + + Returns: + :class:`ETESMILPResult` containing optimal sizes, dispatch + time-series, and cost breakdown. + """ + price = np.asarray(electricity_price, dtype=float) + load = np.asarray(heat_load_kW, dtype=float) + if price.shape != load.shape or price.ndim != 1: + raise ValueError( + f"electricity_price and heat_load_kW must be 1-D arrays of the " + f"same length; got {price.shape} and {load.shape}." + ) + n = len(load) + T = range(n) + + cfg = config + + m = pyo.ConcreteModel() + m.T = pyo.Set(initialize=list(T), ordered=True) + + # --- Sizing variables --- + m.S_TES = pyo.Var( + domain=pyo.NonNegativeReals, + bounds=(cfg.S_TES_min_kWh, cfg.S_TES_max_kWh), + initialize=cfg.S_TES_fixed_kWh + if cfg.S_TES_fixed_kWh is not None + else 0.5 * (cfg.S_TES_min_kWh + cfg.S_TES_max_kWh), + ) + if cfg.S_TES_fixed_kWh is not None: + m.S_TES.fix(cfg.S_TES_fixed_kWh) + + if cfg.etes_type == "P-ETES": + m.S_ch = pyo.Var( + domain=pyo.NonNegativeReals, + bounds=(cfg.S_ch_min_kW, cfg.S_ch_max_kW), + initialize=cfg.S_ch_fixed_kW + if cfg.S_ch_fixed_kW is not None + else 0.5 * (cfg.S_ch_min_kW + cfg.S_ch_max_kW), + ) + if cfg.S_ch_fixed_kW is not None: + m.S_ch.fix(cfg.S_ch_fixed_kW) + m.S_dis = pyo.Var( + domain=pyo.NonNegativeReals, + bounds=(cfg.S_dis_min_kW, cfg.S_dis_max_kW), + initialize=cfg.S_dis_fixed_kW + if cfg.S_dis_fixed_kW is not None + else 0.5 * (cfg.S_dis_min_kW + cfg.S_dis_max_kW), + ) + if cfg.S_dis_fixed_kW is not None: + m.S_dis.fix(cfg.S_dis_fixed_kW) + else: + # R-ETES: S_ch and S_dis are not sized; fix to zero (cost = 0) + m.S_ch = pyo.Var(domain=pyo.NonNegativeReals, initialize=0.0) + m.S_ch.fix(0.0) + m.S_dis = pyo.Var(domain=pyo.NonNegativeReals, initialize=0.0) + m.S_dis.fix(0.0) + + # --- Dispatch variables --- + m.E_st = pyo.Var(m.T, domain=pyo.NonNegativeReals) + m.Q_ch = pyo.Var(m.T, domain=pyo.NonNegativeReals) + m.Q_dis = pyo.Var(m.T, domain=pyo.NonNegativeReals) + m.unmet = pyo.Var(m.T, domain=pyo.NonNegativeReals) + if not cfg.allow_unmet_load: + for t in m.T: + m.unmet[t].fix(0.0) + + # Optional binary on/off vars for mutual exclusion + if cfg.enforce_mutex_charge_dis: + m.z_ch = pyo.Var(m.T, domain=pyo.Binary) + m.z_dis = pyo.Var(m.T, domain=pyo.Binary) + # Big-M based on storage capacity / rate caps + bigM = max(cfg.S_TES_max_kWh / max(dt_h, 1e-9), cfg.S_ch_max_kW, cfg.S_dis_max_kW) + + def _ch_on_rule(mm, t): + return mm.Q_ch[t] <= bigM * mm.z_ch[t] + + def _dis_on_rule(mm, t): + return mm.Q_dis[t] <= bigM * mm.z_dis[t] + + def _mutex_rule(mm, t): + return mm.z_ch[t] + mm.z_dis[t] <= 1 + + m.ch_on = pyo.Constraint(m.T, rule=_ch_on_rule) + m.dis_on = pyo.Constraint(m.T, rule=_dis_on_rule) + m.mutex = pyo.Constraint(m.T, rule=_mutex_rule) + + # --- Storage energy balance --- + E0 = cfg.SOC_init # multiplied by S_TES below + f_loss = cfg.f_loss + + def _balance_rule(mm, t): + # Storage loss energy over the timestep = E_prev * f_loss + # E_st[t] = E_prev*(1 - f_loss) + (Q_ch - Q_dis)*dt + if t == 0: + E_prev = E0 * mm.S_TES + else: + E_prev = mm.E_st[t - 1] + return mm.E_st[t] == E_prev * (1.0 - f_loss) + (mm.Q_ch[t] - mm.Q_dis[t]) * dt_h + + m.balance = pyo.Constraint(m.T, rule=_balance_rule) + + # SOC bounds + def _soc_min_rule(mm, t): + return mm.E_st[t] >= cfg.SOC_min * mm.S_TES + + def _soc_max_rule(mm, t): + return mm.E_st[t] <= mm.S_TES + + m.soc_min = pyo.Constraint(m.T, rule=_soc_min_rule) + m.soc_max = pyo.Constraint(m.T, rule=_soc_max_rule) + + # Cyclic boundary + if cfg.cyclic: + m.cyclic_con = pyo.Constraint(expr=m.E_st[n - 1] == E0 * m.S_TES) + + # Load-meeting constraint: Q_dis * eta_dis + unmet = load + def _load_rule(mm, t): + return mm.Q_dis[t] * cfg.eta_dis + mm.unmet[t] == load[t] + + m.load_con = pyo.Constraint(m.T, rule=_load_rule) + + # Rate / sizing limits per ETES type + if cfg.etes_type == "P-ETES": + # Q_ch <= S_ch * eta_ch (charging unit cap; Q_ch is thermal added) + def _ch_cap_rule(mm, t): + return mm.Q_ch[t] <= mm.S_ch * cfg.eta_ch + + # Q_dis <= S_dis / eta_dis (discharging unit cap; thermal output S_dis) + def _dis_cap_rule(mm, t): + return mm.Q_dis[t] * cfg.eta_dis <= mm.S_dis + + m.ch_cap = pyo.Constraint(m.T, rule=_ch_cap_rule) + m.dis_cap = pyo.Constraint(m.T, rule=_dis_cap_rule) + + # Minimum charging time: S_ch * eta_ch <= S_TES / t_ch_min + if cfg.t_ch_min_h > 0: + m.t_ch_min_con = pyo.Constraint(expr=m.S_ch * cfg.eta_ch * cfg.t_ch_min_h <= m.S_TES) + else: # R-ETES: rates coupled to S_TES + + def _ch_cap_rule(mm, t): + return mm.Q_ch[t] <= cfg.f_ch_max * mm.S_TES + + def _dis_cap_rule(mm, t): + return mm.Q_dis[t] <= cfg.f_dis_max * mm.S_TES + + m.ch_cap = pyo.Constraint(m.T, rule=_ch_cap_rule) + m.dis_cap = pyo.Constraint(m.T, rule=_dis_cap_rule) + + # --- Cost variables --- + # Storage cost: c_tes >= linear; c_tes >= floor (handles economy-of-scale floor) + m.c_tes = pyo.Var(domain=pyo.NonNegativeReals) + m.c_tes_lin = pyo.Constraint(expr=m.c_tes >= cfg.C_lin_TES * m.S_TES + cfg.C_const_TES) + if cfg.C_min_TES > 0: + m.c_tes_floor = pyo.Constraint(expr=m.c_tes >= cfg.C_min_TES * m.S_TES) + + # Charging / discharging unit costs (P-ETES only; zero for R-ETES since S_ch=S_dis=0) + m.c_ch = pyo.Var(domain=pyo.NonNegativeReals) + m.c_dis = pyo.Var(domain=pyo.NonNegativeReals) + if cfg.etes_type == "P-ETES": + m.c_ch_lin = pyo.Constraint(expr=m.c_ch >= cfg.C_lin_ch * m.S_ch + cfg.C_const_ch) + m.c_dis_lin = pyo.Constraint(expr=m.c_dis >= cfg.C_lin_dis * m.S_dis + cfg.C_const_dis) + else: + m.c_ch.fix(0.0) + m.c_dis.fix(0.0) + + # --- Objective: total annualized cost --- + # CapEx = c_tes + c_ch + c_dis + # Annualized capex = (FCR + opex_fraction) * CapEx + # Electricity cost = sum_t (Q_ch[t] / eta_ch) * price[t] * dt + annualization = cfg.fixed_charge_rate + cfg.opex_fraction + grid_cost = sum((m.Q_ch[t] / cfg.eta_ch) * float(price[t]) * dt_h for t in m.T) + unmet_penalty = sum(m.unmet[t] * cfg.unmet_load_penalty * dt_h for t in m.T) + + m.obj = pyo.Objective( + expr=annualization * (m.c_tes + m.c_ch + m.c_dis) + grid_cost + unmet_penalty, + sense=pyo.minimize, + ) + + # --- Solve --- + solver, solver_name = _pick_solver(cfg.solver) + if cfg.solver_options: + for k, v in cfg.solver_options.items(): + solver.options[k] = v + results = solver.solve(m, tee=False) + status = str(results.solver.status) + term = str(results.solver.termination_condition) + if term.lower() not in ("optimal", "feasible", "globally_optimal", "locallyoptimal"): + raise RuntimeError( + f"Solver {solver_name!r} did not find an optimal solution: " + f"status={status}, termination={term}" + ) + + # --- Extract results --- + S_TES = float(pyo.value(m.S_TES)) + S_ch = float(pyo.value(m.S_ch)) + S_dis = float(pyo.value(m.S_dis)) + + E_st = np.array([pyo.value(m.E_st[t]) for t in T]) + Q_ch = np.array([pyo.value(m.Q_ch[t]) for t in T]) + Q_dis = np.array([pyo.value(m.Q_dis[t]) for t in T]) + unmet = np.array([pyo.value(m.unmet[t]) for t in T]) + + # Reconstruct loss series consistent with the formulation + E_prev_series = np.empty(n) + E_prev_series[0] = cfg.SOC_init * S_TES + E_prev_series[1:] = E_st[:-1] + Q_st_loss = E_prev_series * cfg.f_loss / dt_h if dt_h > 0 else np.zeros(n) + Q_ch_loss = Q_ch * (1.0 - cfg.eta_ch) / cfg.eta_ch if cfg.eta_ch > 0 else np.zeros(n) + Q_dis_loss = Q_dis * (1.0 - cfg.eta_dis) + P_grid = Q_ch / cfg.eta_ch if cfg.eta_ch > 0 else np.zeros(n) + + capex = float(pyo.value(m.c_tes + m.c_ch + m.c_dis)) + annual_opex = cfg.opex_fraction * capex + annualized_capex = cfg.fixed_charge_rate * capex + annual_electricity_cost = float(np.sum(P_grid * price * dt_h)) + total = annualized_capex + annual_opex + annual_electricity_cost + + return ETESMILPResult( + S_TES_kWh=S_TES, + S_ch_kW=S_ch, + S_dis_kW=S_dis, + E_st_kWh=E_st, + Q_ch_kW=Q_ch, + Q_dis_kW=Q_dis, + Q_st_loss_kW=Q_st_loss, + Q_ch_loss_kW=Q_ch_loss, + Q_dis_loss_kW=Q_dis_loss, + P_grid_kW=P_grid, + unmet_load_kW=unmet, + capex_USD=capex, + annual_opex_USD=annual_opex, + annual_electricity_cost_USD=annual_electricity_cost, + annualized_capex_USD=annualized_capex, + total_annualized_cost_USD=total, + solver_status=status, + termination_condition=term, + objective_value=float(pyo.value(m.obj)), + ) diff --git a/h2integrate/converters/heat/heat_exchanger_model/hx_shell_tube_steady.py b/h2integrate/converters/heat/heat_exchanger_model/hx_shell_tube_steady.py new file mode 100644 index 000000000..59abfba1e --- /dev/null +++ b/h2integrate/converters/heat/heat_exchanger_model/hx_shell_tube_steady.py @@ -0,0 +1,916 @@ +""" +Created on Wed Dec 10 09:34:19 2025 + +@author: psaini +""" + +from math import log + +import numpy as np +from CoolProp.CoolProp import PropsSI + + +# ---------------------------------------------------------------------- +# 1) CoolProp-based generic fluid properties +# ---------------------------------------------------------------------- +def make_coolprop_fluid(fluid_name="Water", P_bar=1.0): + """ + Returns a function fluid_fun(T_c) that uses CoolProp for fluid properties. + + Parameters + ---------- + fluid_name : str + CoolProp fluid name, e.g. 'Water', 'Ethanol', 'R134a', 'MEG', etc. + P_bar : float + Operating pressure in bar (absolute). Converted internally to Pa. + + Returns + ------- + fluid_fun : callable + fluid_fun(T_c) -> dict(rho, mu, cp, k, Pr) + where T_c is in degC and properties are SI units: + - rho : density [kg/m^3] + - mu : dynamic viscosity [Pa·s] + - cp : specific heat at constant pressure [J/kg-K] + - k : thermal conductivity [W/m-K] + - Pr : Prandtl number [-] + """ + P_Pa = P_bar * 1e5 # bar -> Pa + + def fluid_fun(T_c): + T_c_arr = np.asarray(T_c, dtype=float) + T_K = T_c_arr + 273.15 + + rho = PropsSI("D", "T", T_K, "P", P_Pa, fluid_name) # density [kg/m^3] + mu = PropsSI("V", "T", T_K, "P", P_Pa, fluid_name) # dyn. visc. [Pa·s] + cp_ = PropsSI("C", "T", T_K, "P", P_Pa, fluid_name) # cp [J/kg-K] + k = PropsSI("L", "T", T_K, "P", P_Pa, fluid_name) # k [W/m-K] + Pr = PropsSI("Prandtl", "T", T_K, "P", P_Pa, fluid_name) # Pr [-] + + return {"rho": rho, "mu": mu, "cp": cp_, "k": k, "Pr": Pr} + + return fluid_fun + + +# Default fluids (both sides = liquid water @ 1 bar) +DEFAULT_WATER_FLUID = make_coolprop_fluid("Water", P_bar=1.0) + + +# ---------------------------------------------------------------------- +# 2) Metal properties, helper functions +# ---------------------------------------------------------------------- +def metal_props(name: str): + key = name.lower().replace(" ", "_") + if key in ("carbon_steel", "cs"): + return 45.0, 7850.0, 480.0 + elif key in ("stainless_304", "ss304", "stainless"): + return 16.0, 8000.0, 500.0 + elif key in ("copper", "cu"): + return 390.0, 8960.0, 385.0 + else: + # default carbon steel + return 45.0, 7850.0, 480.0 + + +def LMTD_calc(dT1, dT2): + if dT1 <= 0 or dT2 <= 0: + raise ValueError("Non-positive temperature difference in LMTD_calc.") + if abs(dT1 - dT2) < 1e-6: + return 0.5 * (dT1 + dT2) + return (dT1 - dT2) / log(dT1 / dT2) + + +def classify_regime(Re_min, Re_max, Re_laminar_max, Re_turb_min): + if Re_max < Re_laminar_max: + return "laminar" + elif Re_min > Re_turb_min: + return "turbulent" + else: + return "transition" + + +def shell_h_local(cold_props, m_dot_c, geom, shell, F_shell, model): + S_m = geom["S_m"] + D_o = geom["D_o"] + p_t = shell["pitch"] + D_e = 4.0 * (p_t**2 - np.pi * D_o**2 / 4.0) / (np.pi * D_o) + + G_s = m_dot_c / S_m + Re_s = G_s * D_e / cold_props["mu"] + Pr_c = cold_props["cp"] * cold_props["mu"] / cold_props["k"] + + model_l = model.lower() + if model_l == "kern": + if Re_s < 2000.0: + L_char = max(1e-6, shell["baffle_spacing"]) + Nu_fd = 3.66 + Gz_s = Re_s * Pr_c * D_e / L_char + Nu_s = Nu_fd * (1.0 - np.exp(-Gz_s / 20.0)) + Nu_s = max(1.0, min(Nu_fd, Nu_s)) + else: + Nu_s = 0.36 * Re_s**0.55 * Pr_c ** (1.0 / 3.0) + elif model_l == "bell_delaware": + # Placeholder: same as Kern, can refine later + if Re_s < 2000.0: + L_char = max(1e-6, shell["baffle_spacing"]) + Nu_fd = 3.66 + Gz_s = Re_s * Pr_c * D_e / L_char + Nu_s = Nu_fd * (1.0 - np.exp(-Gz_s / 20.0)) + Nu_s = max(1.0, min(Nu_fd, Nu_s)) + else: + Nu_s = 0.36 * Re_s**0.55 * Pr_c ** (1.0 / 3.0) + else: + raise ValueError(f"Unknown shell_model '{model}'") + + h_o = Nu_s * cold_props["k"] / D_e + h_o *= F_shell + return h_o, Re_s, Nu_s + + +def shell_dp_Kern(m_dot_c, cold_props, geom, shell): + S_m = geom["S_m"] + D_o = geom["D_o"] + p_t = shell["pitch"] + D_e = 4.0 * (p_t**2 - np.pi * D_o**2 / 4.0) / (np.pi * D_o) + + G_s = m_dot_c / S_m + Re_s = G_s * D_e / cold_props["mu"] + + if Re_s < 2000.0: + f_s = 64.0 / Re_s + else: + f_s = 0.14 * Re_s ** (-0.2) + + N_b = max(1, int(geom["L_tube"] / shell["baffle_spacing"]) - 1) + dp_shell = 4.0 * f_s * (G_s**2 / (2.0 * cold_props["rho"])) * N_b + return dp_shell + + +# ---------------------------------------------------------------------- +# 3) Main HX model +# ---------------------------------------------------------------------- +def hx_shell_tube_steady(params=None): + if params is None: + params = {} + get = params.get + + # Inlet temps [°C] and mass flow rates [kg/s] + Th_in = float(get("Th_in", 90.0)) + Tc_in = float(get("Tc_in", 30.0)) + m_dot_h = float(get("m_dot_h", 18.0)) + m_dot_c = float(get("m_dot_c", 40.0)) + + # Geometry + N_tubes = int(get("N_tubes", 192)) + N_passes = int(get("N_passes", 2)) + L_tube = float(get("L_tube", 4.9)) + D_o = float(get("D_o", 0.01905)) + t_wall = float(get("t_wall", 0.002)) + D_i = D_o - 2.0 * t_wall + k_wall = float(get("k_wall", 45.0)) + D_shell = float(get("D_shell", 0.591)) + + # Discretization + N_seg = int(get("N_seg", 20)) + + # Fouling + R_fi = float(get("R_fi", 1e-4)) + R_fo = float(get("R_fo", 2e-4)) + + # Shell-side parameters + shell = { + "pitch": get("pitch", 1.1 * D_o), + "baffle_spacing": get("baffle_spacing", 0.20), + "baffle_cut": get("baffle_cut", 0.20), + "layout": get("layout", "triangular"), + } + + # Pump efficiency + eta_pump = float(get("eta_pump", 0.7)) + + # Pinch threshold + pinch_threshold = float(get("pinch_threshold", 3.0)) + + # Physics switches + use_SiederTate_tube = bool(get("use_SiederTate_tube", True)) + use_wall_viscosity_correction = bool(get("use_wall_viscosity_correction", True)) + F_shell = float(get("F_shell", 0.8)) + shell_model = get("shell_model", "kern") + T0_env = float(get("T0_env", 298.15)) + + # External shell loss + U_loss = float(get("U_loss", 0.0)) + T_amb = float(get("T_amb", 25.0)) + + # Property functions (now CoolProp-based by default) + fluid_hot_fun = get("fluid_hot", DEFAULT_WATER_FLUID) + fluid_cold_fun = get("fluid_cold", DEFAULT_WATER_FLUID) + + # fzero options + fzero_tol = float(get("fzero_tol", 1e-6)) + fzero_maxit = int(get("fzero_maxit", 50)) + + # Materials + tube_material = get("tube_material", "carbon_steel") + shell_material = get("shell_material", "carbon_steel") + + t_shell = float(get("t_shell", 0.02)) + + # Reynolds thresholds + Re_laminar_max = float(get("Re_laminar_max", 2300.0)) + Re_turb_min = float(get("Re_turb_min", 4000.0)) + + # Metal props + tube_k, tube_rho, tube_cp = metal_props(tube_material) + shell_k, shell_rho, shell_cp = metal_props(shell_material) + + if "k_wall" not in params: + k_wall = tube_k + + # Geometry and areas + A_tube_flow = (np.pi * D_i**2 / 4.0) * (N_tubes / N_passes) + A_o_total = np.pi * D_o * L_tube * N_tubes + A_seg = A_o_total / N_seg + + A_shell_free = (np.pi * D_shell**2 / 4.0) - N_tubes * (np.pi * D_o**2 / 4.0) + if A_shell_free <= 0.0: + raise ValueError("Shell is too small for the given number/size of tubes.") + + P_wet_shell = np.pi * D_shell + N_tubes * np.pi * D_o + D_h_shell = 4.0 * A_shell_free / P_wet_shell + + D_shell_in = D_shell + D_shell_out = D_shell + 2.0 * t_shell + + V_tube_metal = (np.pi / 4.0) * (D_o**2 - D_i**2) * L_tube * N_tubes + V_shell_metal = (np.pi / 4.0) * (D_shell_out**2 - D_shell_in**2) * L_tube + + V_tube_fluid = (np.pi / 4.0) * D_i**2 * L_tube * N_tubes + V_shell_fluid = A_shell_free * L_tube + + S_m = D_shell * shell["baffle_spacing"] * (1.0 - shell["baffle_cut"]) + + geom = { + "N_tubes": N_tubes, + "N_passes": N_passes, + "L_tube": L_tube, + "D_i": D_i, + "D_o": D_o, + "D_shell": D_shell, + "A_tube_flow": A_tube_flow, + "A_shell_free": A_shell_free, + "D_h_shell": D_h_shell, + "A_o_total": A_o_total, + "A_seg": A_seg, + "N_seg": N_seg, + "A_shell_outer_total": np.pi * D_shell * L_tube, + "S_m": S_m, + } + geom["A_shell_outer_seg"] = geom["A_shell_outer_total"] / N_seg + + # ------------------------------------------------------------------ + # Shooting for counter-current + # ------------------------------------------------------------------ + Tc_out0_min = Tc_in + 1e-3 + Tc_out0_max = Th_in - 1e-3 + Tc_out0_guess = Tc_in + 0.7 * (Th_in - Tc_in) + + def march_segments(Th_in_local, Tc_out0_local): + N_seg = geom["N_seg"] + A_seg = geom["A_seg"] + A_tube_flow = geom["A_tube_flow"] + D_i = geom["D_i"] + + Th = np.zeros(N_seg + 1) + Tc = np.zeros(N_seg + 1) + Tw = np.zeros(N_seg) + U_seg = np.zeros(N_seg) + Q_seg = np.zeros(N_seg) + Q_loss_seg = np.zeros(N_seg) + Re_t_seg = np.zeros(N_seg) + Re_s_seg = np.zeros(N_seg) + h_i_seg = np.zeros(N_seg) + h_o_seg = np.zeros(N_seg) + Nu_t_seg = np.zeros(N_seg) + Nu_s_seg = np.zeros(N_seg) + S_gen_seg = np.zeros(N_seg) + + Th[0] = Th_in_local + Tc[0] = Tc_out0_local + + A_loss = geom["A_shell_outer_seg"] + L_seg = geom["L_tube"] / geom["N_seg"] + + for i in range(N_seg): + Th_i = Th[i] + Tc_i = Tc[i] + + hot = fluid_hot_fun(Th_i) + cold = fluid_cold_fun(Tc_i) + + # Tube-side Re, Nu, h_i + V_tube = m_dot_h / (hot["rho"] * A_tube_flow) + Re_t = hot["rho"] * V_tube * D_i / hot["mu"] + Pr_h = hot["cp"] * hot["mu"] / hot["k"] + + if Re_t < Re_laminar_max: + Nu_fd = 3.66 + Gz_t = Re_t * Pr_h * D_i / max(1e-9, L_seg) + Nu_t = Nu_fd * (1.0 - np.exp(-Gz_t / 20.0)) + Nu_t = max(1.0, min(Nu_fd, Nu_t)) + else: + if use_SiederTate_tube: + mu_b = hot["mu"] + if use_wall_viscosity_correction: + Tw_guess = 0.5 * (Th_i + Tc_i) + wall_hot = fluid_hot_fun(Tw_guess) + mu_w = wall_hot["mu"] + else: + mu_w = mu_b + Nu_t = 0.027 * Re_t**0.8 * Pr_h ** (1.0 / 3.0) * (mu_b / mu_w) ** 0.14 + else: + Nu_t = 0.023 * Re_t**0.8 * Pr_h**0.4 + h_i = Nu_t * hot["k"] / D_i + + Re_t_seg[i] = Re_t + h_i_seg[i] = h_i + Nu_t_seg[i] = Nu_t + + # Shell-side + h_o, Re_s, Nu_s = shell_h_local(cold, m_dot_c, geom, shell, F_shell, shell_model) + Re_s_seg[i] = Re_s + h_o_seg[i] = h_o + Nu_s_seg[i] = Nu_s + + # Overall U + R_o = 1.0 / h_o + R_fo + R_w = (geom["D_o"] * log(geom["D_o"] / geom["D_i"])) / (2.0 * k_wall) + R_i = geom["D_o"] / (geom["D_i"] * h_i) + R_fi * geom["D_o"] / geom["D_i"] + + U_i = 1.0 / (R_o + R_w + R_i) + U_seg[i] = U_i + + dT_i = Th_i - Tc_i + Q_i = U_i * A_seg * dT_i + + # External loss (from shell to ambient) + Q_loss_i = U_loss * A_loss * (Tc_i - T_amb) + Q_loss_seg[i] = Q_loss_i + + cp_h = hot["cp"] + cp_c = cold["cp"] + + # Energy balances + Th[i + 1] = Th_i - Q_i / (m_dot_h * cp_h) + Tc[i + 1] = Tc_i - (Q_i - Q_loss_i) / (m_dot_c * cp_c) + + # Wall temperature + Tw[i] = (h_i * Th_i + h_o * Tc_i) / (h_i + h_o) + Q_seg[i] = Q_i + + # Local entropy generation + Th_iK = Th_i + 273.15 + Tc_iK = Tc_i + 273.15 + S_gen_seg[i] = Q_i * (1.0 / Tc_iK - 1.0 / Th_iK) + + return ( + Th, + Tc, + Tw, + U_seg, + Q_seg, + Q_loss_seg, + Re_t_seg, + Re_s_seg, + h_i_seg, + h_o_seg, + Nu_t_seg, + Nu_s_seg, + S_gen_seg, + ) + + def residual_Tc_out(Tc_out0): + Th_, Tc_, *_ = march_segments(Th_in, Tc_out0) + Tc_at_L = Tc_[-1] + return Tc_at_L - Tc_in + + # Simple bisection + def find_root_bisect(f, a, b, tol, maxit): + fa = f(a) + fb = f(b) + if fa * fb > 0: + return None, fa, fb + for _ in range(maxit): + c = 0.5 * (a + b) + fc = f(c) + if abs(fc) < tol or 0.5 * (b - a) < tol: + return c, fa, fb + if fa * fc < 0: + b = c + fb = fc + else: + a = c + fa = fc + return c, fa, fb + + Tc_root, r_min, r_max = find_root_bisect( + residual_Tc_out, Tc_out0_min, Tc_out0_max, fzero_tol, fzero_maxit + ) + shooting_bracket_ok = True + if Tc_root is None: + shooting_bracket_ok = False + Tc_root = Tc_out0_guess + + # March with final + ( + Th, + Tc, + Tw, + U_seg, + Q_seg, + Q_loss_seg, + Re_t_seg, + Re_s_seg, + h_i_seg, + h_o_seg, + Nu_t_seg, + Nu_s_seg, + S_gen_seg, + ) = march_segments(Th_in, Tc_root) + + x = np.linspace(0.0, L_tube, N_seg + 1) + x_mid = 0.5 * (x[:-1] + x[1:]) + + Q_total = float(np.sum(Q_seg)) + Q_loss_tot = float(np.sum(Q_loss_seg)) + Th_out = float(Th[-1]) + Tc_out = float(Tc[0]) + + # Global metrics + dT1 = Th_in - Tc_out + dT2 = Th_out - Tc_in + LMTD_global = LMTD_calc(dT1, dT2) + U_global = Q_total / (A_o_total * LMTD_global) + + hot_in = fluid_hot_fun(Th_in) + cold_in = fluid_cold_fun(Tc_in) + cp_h_in = hot_in["cp"] + cp_c_in = cold_in["cp"] + C_hot = m_dot_h * cp_h_in + C_cold = m_dot_c * cp_c_in + C_min = min(C_hot, C_cold) + C_max = max(C_hot, C_cold) + C_r = C_min / C_max + + epsilon = Q_total / (C_min * (Th_in - Tc_in)) + NTU = U_global * A_o_total / C_min + + # Energy-balance checks + Q_hot = m_dot_h * cp_h_in * (Th_in - Th_out) + Q_cold = m_dot_c * cp_c_in * (Tc_out - Tc_in) + + Q_ref_hot = Q_total + Q_ref_cold = Q_total - Q_loss_tot + + energy_balance_error_hot = (Q_hot - Q_ref_hot) / abs(Q_ref_hot) + energy_balance_error_cold = (Q_cold - Q_ref_cold) / abs(Q_ref_cold) + + # Lumped U & LMTD check + Th_mean_bulk = 0.5 * (Th_in + Th_out) + Tc_mean_bulk = 0.5 * (Tc_in + Tc_out) + hot_m = fluid_hot_fun(Th_mean_bulk) + cold_m = fluid_cold_fun(Tc_mean_bulk) + Tw_mean = 0.5 * (Th_mean_bulk + Tc_mean_bulk) + + V_tube_mean = m_dot_h / (hot_m["rho"] * A_tube_flow) + Re_t_mean = hot_m["rho"] * V_tube_mean * D_i / hot_m["mu"] + Pr_h_mean = hot_m["cp"] * hot_m["mu"] / hot_m["k"] + + if Re_t_mean < Re_laminar_max: + Nu_t_mean = 3.66 + else: + if use_SiederTate_tube: + mu_b = hot_m["mu"] + if use_wall_viscosity_correction: + wall_hot_mean = fluid_hot_fun(Tw_mean) + mu_w = wall_hot_mean["mu"] + else: + mu_w = mu_b + Nu_t_mean = 0.027 * Re_t_mean**0.8 * Pr_h_mean ** (1.0 / 3.0) * (mu_b / mu_w) ** 0.14 + else: + Nu_t_mean = 0.023 * Re_t_mean**0.8 * Pr_h_mean**0.4 + h_i_mean = Nu_t_mean * hot_m["k"] / D_i + + h_o_mean, _, _ = shell_h_local(cold_m, m_dot_c, geom, shell, F_shell, shell_model) + + R_o_mean = 1.0 / h_o_mean + R_fo + R_w_mean = (D_o * log(D_o / D_i)) / (2.0 * k_wall) + R_i_mean = D_o / (D_i * h_i_mean) + R_fi * D_o / D_i + + U_lumped = 1.0 / (R_o_mean + R_w_mean + R_i_mean) + Q_LMTD_lumped = U_lumped * A_o_total * LMTD_global + Q_error_frac = (Q_total - Q_LMTD_lumped) / Q_total + + # Pressure drops & pump power + Th_mean = float(np.mean(Th)) + Tc_mean = float(np.mean(Tc)) + hot_mean = fluid_hot_fun(Th_mean) + cold_mean = fluid_cold_fun(Tc_mean) + + V_tube = m_dot_h / (hot_mean["rho"] * A_tube_flow) + Re_t = hot_mean["rho"] * V_tube * D_i / hot_mean["mu"] + + if Re_t < Re_laminar_max: + f_t = 64.0 / Re_t + else: + f_t = 0.3164 * Re_t ** (-0.25) + + L_equiv_tube = L_tube * N_passes + dp_fric_tube = f_t * (L_equiv_tube / D_i) * 0.5 * hot_mean["rho"] * V_tube**2 + + K_inlet = 1.0 + K_exit = 1.0 + K_bends = 1.5 * (N_passes - 1) + K_total = K_inlet + K_exit + K_bends + dp_minor_tube = K_total * 0.5 * hot_mean["rho"] * V_tube**2 + + dp_tube_total = dp_fric_tube + dp_minor_tube + + dp_shell_total = shell_dp_Kern(m_dot_c, cold_mean, geom, shell) + + P_pump_tube = (m_dot_h / hot_mean["rho"]) * dp_tube_total / eta_pump + P_pump_shell = (m_dot_c / cold_mean["rho"]) * dp_shell_total / eta_pump + P_pump_total = P_pump_tube + P_pump_shell + + # Inventories + time constants + rho_hot_mean = hot_mean["rho"] + rho_cold_mean = cold_mean["rho"] + cp_hot_mean = hot_mean["cp"] + cp_cold_mean = cold_mean["cp"] + + m_tube_metal = tube_rho * V_tube_metal + m_shell_metal = shell_rho * V_shell_metal + m_tube_fluid = rho_hot_mean * V_tube_fluid + m_shell_fluid = rho_cold_mean * V_shell_fluid + + C_tube_metal = m_tube_metal * tube_cp + C_shell_metal = m_shell_metal * shell_cp + C_tube_fluid = m_tube_fluid * cp_hot_mean + C_shell_fluid = m_shell_fluid * cp_cold_mean + + C_hot_total = C_tube_metal + C_tube_fluid + C_cold_total = C_shell_metal + C_shell_fluid + + tau_hot = C_hot_total / (U_global * A_o_total) + tau_cold = C_cold_total / (U_global * A_o_total) + + # Pinch & crossing + deltaT = Th - Tc + min_deltaT = float(np.min(deltaT)) + has_cross = bool(np.any(deltaT <= 0.0)) + pinch_warning = (not has_cross) and (min_deltaT < pinch_threshold) + + # 2nd law metrics + Th_in_K = Th_in + 273.15 + Th_out_K = Th_out + 273.15 + Tc_in_K = Tc_in + 273.15 + Tc_out_K = Tc_out + 273.15 + + hot_in2 = fluid_hot_fun(Th_in) + hot_out2 = fluid_hot_fun(Th_out) + cold_in2 = fluid_cold_fun(Tc_in) + cold_out2 = fluid_cold_fun(Tc_out) + + cp_h_bar = 0.5 * (hot_in2["cp"] + hot_out2["cp"]) + cp_c_bar = 0.5 * (cold_in2["cp"] + cold_out2["cp"]) + + S_gen_dot_bulk = m_dot_h * cp_h_bar * log(Th_out_K / Th_in_K) + m_dot_c * cp_c_bar * log( + Tc_out_K / Tc_in_K + ) + Ex_dest_dot_bulk = T0_env * S_gen_dot_bulk + + S_gen_dot_seg = float(np.sum(S_gen_seg)) + Ex_dest_dot_seg = T0_env * S_gen_dot_seg + + # Regimes + Re_tube_min = float(np.min(Re_t_seg)) + Re_tube_max = float(np.max(Re_t_seg)) + Re_shell_min = float(np.min(Re_s_seg)) + Re_shell_max = float(np.max(Re_s_seg)) + + tube_regime = classify_regime(Re_tube_min, Re_tube_max, Re_laminar_max, Re_turb_min) + shell_regime = classify_regime(Re_shell_min, Re_shell_max, Re_laminar_max, Re_turb_min) + + V_tube_mean2 = V_tube + G_s_mean = m_dot_c / S_m + V_shell_mean = G_s_mean / cold_mean["rho"] + + flags = { + "has_cross": has_cross, + "small_pinch": pinch_warning, + "tube_laminar": (tube_regime == "laminar"), + "shell_laminar": (shell_regime == "laminar"), + "low_NTU": (NTU < 0.1), + "high_dp_tube": (dp_tube_total > 1e5), + "high_dp_shell": (dp_shell_total > 1e4), + "Q_LMTD_mismatch": (abs(Q_error_frac) > 0.05), + "shooting_ok": shooting_bracket_ok, + "V_tube_low": (V_tube_mean2 < 1.0), + "V_tube_high": (V_tube_mean2 > 3.0), + "V_shell_low": (V_shell_mean < 0.3), + "V_shell_high": (V_shell_mean > 1.5), + } + + inventory = { + "tube_material": tube_material, + "shell_material": shell_material, + "t_shell": t_shell, + "V_tube_metal": V_tube_metal, + "V_shell_metal": V_shell_metal, + "V_tube_fluid": V_tube_fluid, + "V_shell_fluid": V_shell_fluid, + "m_tube_metal": m_tube_metal, + "m_shell_metal": m_shell_metal, + "m_tube_fluid": m_tube_fluid, + "m_shell_fluid": m_shell_fluid, + "C_tube_metal": C_tube_metal, + "C_shell_metal": C_shell_metal, + "C_tube_fluid": C_tube_fluid, + "C_shell_fluid": C_shell_fluid, + "C_hot_total": C_hot_total, + "C_cold_total": C_cold_total, + } + + res = { + "x": x, + "x_mid": x_mid, + "Th": Th, + "Tc": Tc, + "Tw": Tw, + "Q_seg": Q_seg, + "Q_total": Q_total, + "Q_loss_seg": Q_loss_seg, + "Q_loss_tot": Q_loss_tot, + "U_seg": U_seg, + "U_global": U_global, + "LMTD_global": LMTD_global, + "epsilon": epsilon, + "NTU": NTU, + "C_r": C_r, + "Q_hot": Q_hot, + "Q_cold": Q_cold, + "energy_balance_error_hot": energy_balance_error_hot, + "energy_balance_error_cold": energy_balance_error_cold, + "dp_tube_total": dp_tube_total, + "dp_shell_total": dp_shell_total, + "P_pump_tube": P_pump_tube, + "P_pump_shell": P_pump_shell, + "P_pump_total": P_pump_total, + "min_deltaT": min_deltaT, + "has_cross": has_cross, + "pinch_warning": pinch_warning, + "Re_tube": Re_t_seg, + "Re_shell": Re_s_seg, + "Re_tube_min": Re_tube_min, + "Re_tube_max": Re_tube_max, + "Re_shell_min": Re_shell_min, + "Re_shell_max": Re_shell_max, + "tube_regime": tube_regime, + "shell_regime": shell_regime, + "h_i_seg": h_i_seg, + "h_o_seg": h_o_seg, + "Nu_t_seg": Nu_t_seg, + "Nu_s_seg": Nu_s_seg, + "U_lumped": U_lumped, + "Q_LMTD_lumped": Q_LMTD_lumped, + "Q_error_frac": Q_error_frac, + "S_gen_seg": S_gen_seg, + "S_gen_dot": S_gen_dot_seg, + "Ex_dest_dot": Ex_dest_dot_seg, + "S_gen_dot_bulk": S_gen_dot_bulk, + "Ex_dest_dot_bulk": Ex_dest_dot_bulk, + "tau_hot": tau_hot, + "tau_cold": tau_cold, + "inventory": inventory, + "V_tube_mean": V_tube_mean2, + "V_shell_mean": V_shell_mean, + "flags": flags, + "geom": geom, + "shell": shell, + "params": { + "Th_in": Th_in, + "Tc_in": Tc_in, + "m_dot_h": m_dot_h, + "m_dot_c": m_dot_c, + "N_tubes": N_tubes, + "N_passes": N_passes, + "L_tube": L_tube, + "D_o": D_o, + "t_wall": t_wall, + "D_i": D_i, + "D_shell": D_shell, + "N_seg": N_seg, + "R_fi": R_fi, + "R_fo": R_fo, + "eta_pump": eta_pump, + "pinch_threshold": pinch_threshold, + "use_SiederTate_tube": use_SiederTate_tube, + "use_wall_viscosity_correction": use_wall_viscosity_correction, + "F_shell": F_shell, + "shell_model": shell_model, + "T0_env": T0_env, + "U_loss": U_loss, + "T_amb": T_amb, + "fzero_tol": fzero_tol, + "fzero_maxit": fzero_maxit, + "tube_material": tube_material, + "shell_material": shell_material, + "t_shell": t_shell, + "Re_laminar_max": Re_laminar_max, + "Re_turb_min": Re_turb_min, + }, + } + return res + + +# ---------------------------------------------------------------------- +# 4) Printing summary +# ---------------------------------------------------------------------- +def print_summary(res, label="Python + CoolProp"): + Th_in = res["params"]["Th_in"] + Tc_in = res["params"]["Tc_in"] + Th_out = res["Th"][-1] + Tc_out = res["Tc"][0] + + print(f"=== HX summary (detailed, {label}) ===") + print(f"Hot side: Th_in = {Th_in:6.2f} degC, Th_out = {Th_out:6.2f} degC") + print(f"Cold side: Tc_in = {Tc_in:6.2f} degC, Tc_out = {Tc_out:6.2f} degC\n") + + print(f"Q_total = {res['Q_total']/1e3:8.1f} kW") + print(f"Q_hot (bulk) = {res['Q_hot']/1e3:8.1f} kW") + print(f"Q_cold (bulk) = {res['Q_cold']/1e3:8.1f} kW") + print( + f"Q_loss_tot = {res['Q_loss_tot']/1e3:8.1f} kW " + f"({100*res['Q_loss_tot']/max(1e-9,res['Q_total']):.2f} %)\n" + ) + + print(f"epsilon = {res['epsilon']:.3f}") + print(f"NTU = {res['NTU']:.3f}") + print(f"C_r = {res['C_r']:.3f}") + print(f"LMTD_global = {res['LMTD_global']:.2f} K") + print(f"U_global = {res['U_global']:.1f} W/m^2-K") + print(f"U_lumped = {res['U_lumped']:.1f} W/m^2-K") + print(f"Q_LMTD_lumped = {res['Q_LMTD_lumped']/1e3:8.1f} kW") + print(f"Q_error_frac = {res['Q_error_frac']:.3f} " f"({100*res['Q_error_frac']:.1f} %)\n") + + print(f"U_seg min / max = {res['U_seg'].min():.1f} / {res['U_seg'].max():.1f} W/m^2-K") + print(f"h_i min / max = {res['h_i_seg'].min():.1f} / {res['h_i_seg'].max():.1f} W/m^2-K") + print(f"h_o min / max = {res['h_o_seg'].min():.1f} / {res['h_o_seg'].max():.1f} W/m^2-K\n") + + print( + f"Re_tube min/max = {int(res['Re_tube_min']):6d} / {int(res['Re_tube_max']):6d} " + f"(regime: {res['tube_regime']})" + ) + print( + f"Re_shell min/max = {int(res['Re_shell_min']):6d} / {int(res['Re_shell_max']):6d} " + f"(regime: {res['shell_regime']})" + ) + print(f"V_tube_mean = {res['V_tube_mean']:.3f} m/s") + print(f"V_shell_mean = {res['V_shell_mean']:.3f} m/s\n") + + print(f"dp_tube_total = {res['dp_tube_total']/1e3:.2f} kPa") + print(f"dp_shell_total = {res['dp_shell_total']/1e3:.2f} kPa") + print(f"P_pump_tube = {res['P_pump_tube']/1e3:.3f} kW") + print(f"P_pump_shell = {res['P_pump_shell']/1e3:.3f} kW") + print(f"P_pump_total = {res['P_pump_total']/1e3:.3f} kW\n") + + print(f"tau_hot = {res['tau_hot']:.1f} s") + print(f"tau_cold = {res['tau_cold']:.1f} s\n") + + print(f"min_deltaT = {res['min_deltaT']:.2f} K") + print(f"has_cross = {int(res['has_cross'])}") + print( + f"pinch_warning = {int(res['pinch_warning'])} " + f"(threshold = {res['params']['pinch_threshold']:.1f} K)\n" + ) + + print(f"S_gen_dot (seg) = {res['S_gen_dot']:.2f} W/K") + print(f"Ex_dest_dot (seg)= {res['Ex_dest_dot']/1e3:.2f} kW") + print(f"S_gen_dot_bulk = {res['S_gen_dot_bulk']:.2f} W/K") + print(f"Ex_dest_dot_bulk = {res['Ex_dest_dot_bulk']/1e3:.2f} kW\n") + + print("Flags:") + for k, v in res["flags"].items(): + print(f" {k:15s} = {int(v)}") + print("============================================") + + +# ---------------------------------------------------------------------- +# 5) Plotting +# ---------------------------------------------------------------------- +def plot_results(res): + import matplotlib.pyplot as plt + + x = res["x"] + x_mid = res["x_mid"] + Th = res["Th"] + Tc = res["Tc"] + Tw = res["Tw"] + Q_seg = res["Q_seg"] + S_gen_seg = res["S_gen_seg"] + U_seg = res["U_seg"] + h_i_seg = res["h_i_seg"] + h_o_seg = res["h_o_seg"] + Re_t = res["Re_tube"] + Re_s = res["Re_shell"] + + geom = res["geom"] + L_tube = geom["L_tube"] + N_seg = geom["N_seg"] + L_seg = L_tube / N_seg + + # Heat-transfer per unit length [W/m] + q_per_m = Q_seg / L_seg + + # --- Figure 1: Temperature profiles --- + plt.figure(figsize=(6, 4)) + plt.plot(x, Th, "-r", label="T_h(x)") + plt.plot(x, Tc, "-b", label="T_c(x)") + plt.plot(x_mid, Tw, "--k", label="T_w(x)") + plt.grid(True) + plt.xlabel("Axial position x [m]") + plt.ylabel("Temperature [°C]") + plt.title("Temperature profiles along HX (hot, cold, wall)") + plt.legend() + plt.tight_layout() + + # --- Figure 2: Local q' and S_gen --- + fig, axes = plt.subplots(2, 1, figsize=(6, 6), sharex=True) + axes[0].plot(x_mid, q_per_m / 1e3, "-o") + axes[0].grid(True) + axes[0].set_ylabel("q' [kW/m]") + axes[0].set_title("Local heat transfer per unit length") + + axes[1].plot(x_mid, S_gen_seg, "-s") + axes[1].grid(True) + axes[1].set_xlabel("Axial position x [m]") + axes[1].set_ylabel("S_gen,seg [W/K]") + axes[1].set_title("Local entropy generation per segment") + + fig.tight_layout() + + # --- Figure 3: Local h_i, h_o, and U_seg --- + fig, axes = plt.subplots(2, 1, figsize=(6, 6), sharex=True) + axes[0].plot(x_mid, h_i_seg, "-r", label="h_i (tube side)") + axes[0].plot(x_mid, h_o_seg, "-b", label="h_o (shell side)") + axes[0].grid(True) + axes[0].set_ylabel("h [W/m²-K]") + axes[0].set_title("Local film coefficients") + axes[0].legend() + + axes[1].plot(x_mid, U_seg, "-k") + axes[1].grid(True) + axes[1].set_xlabel("Axial position x [m]") + axes[1].set_ylabel("U [W/m²-K]") + axes[1].set_title("Local overall heat transfer coefficient U(x)") + + fig.tight_layout() + + # --- Figure 4: Reynolds numbers --- + plt.figure(figsize=(6, 4)) + plt.plot(x_mid, Re_t, "-r", label="Re_tube") + plt.plot(x_mid, Re_s, "-b", label="Re_shell") + plt.axhline(res["params"]["Re_laminar_max"], linestyle="--", color="k", label="Re_laminar,max") + plt.axhline(res["params"]["Re_turb_min"], linestyle="--", color="k", label="Re_turb,min") + plt.yscale("log") + plt.grid(True, which="both", axis="y") + plt.xlabel("Axial position x [m]") + plt.ylabel("Re [-]") + plt.title("Local Reynolds numbers (tube & shell)") + plt.legend() + plt.tight_layout() + + +# ---------------------------------------------------------------------- +# 6) Run when executed as a script +# ---------------------------------------------------------------------- +if __name__ == "__main__": + # Example 1: default = Water @ 1 bar on both sides (no need to pass fluids) + # res = hx_shell_tube_steady({}) + + # Example 2: custom CoolProp fluids (e.g., hot water @ 5 bar, cold water @ 2 bar) + hot_fluid = make_coolprop_fluid("Water", P_bar=5.0) + cold_fluid = make_coolprop_fluid("Water", P_bar=2.0) + + params = { + "Th_in": 90.0, + "Tc_in": 30.0, + "m_dot_h": 18.0, + "m_dot_c": 40.0, + "fluid_hot": hot_fluid, + "fluid_cold": cold_fluid, + } + + res = hx_shell_tube_steady(params) + print_summary(res) + plot_results(res) diff --git a/h2integrate/converters/heat/shell_tube_hx.py b/h2integrate/converters/heat/shell_tube_hx.py new file mode 100644 index 000000000..c18507bfa --- /dev/null +++ b/h2integrate/converters/heat/shell_tube_hx.py @@ -0,0 +1,123 @@ +import openmdao.api as om +from attrs import field, define + +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.core.validators import gt_zero +from h2integrate.converters.heat.heat_exchanger_model.hx_shell_tube_steady import ( + hx_shell_tube_steady, +) + + +@define(kw_only=True) +class ShellTubeHXPerformanceModelConfig(BaseConfig): + """ + Configuration class for the ShellTubeHXPerformanceModel. + + Args: + Th_in_C (float): Hot fluid inlet temperature in degrees Celsius. + Tc_in_C (float): Cold fluid inlet temperature in degrees Celsius. + m_dot_h_kg_s (float): Mass flow rate of the hot fluid in kg/s. + m_dot_c_kg_s (float): Mass flow rate of the cold fluid in kg/s. + N_tubes (int): Number of tubes in the heat exchanger. + N_passes (int): Number of passes in the heat exchanger. + L_tube_m (float): Length of each tube in meters. + D_o_m (float): Outer diameter of the tubes in meters. + t_wall_m (float): Wall thickness of the tubes in meters. + D_shell_m (float): Diameter of the shell in meters. + cost_year (int): Year for cost estimation. + """ + + Th_in_C: float = field(validator=gt_zero) + Tc_in_C: float = field(validator=gt_zero) + m_dot_h_kg_s: float = field(validator=gt_zero) + m_dot_c_kg_s: float = field(validator=gt_zero) + N_tubes: int = field(validator=gt_zero) + N_passes: int = field(validator=gt_zero) + L_tube_m: float = field(validator=gt_zero) + D_o_m: float = field(validator=gt_zero) + t_wall_m: float = field(validator=gt_zero) + D_shell_m: float = field(validator=gt_zero) + cost_year: int = field(validator=gt_zero) + + +class ShellTubeHXPerformanceModel(om.ExplicitComponent): + """ + An OpenMDAO component that wraps the steady shell tube heat exchanger model. + """ + + def initialize(self): + self.options.declare("driver_config", types=dict) + self.options.declare("plant_config", types=dict) + self.options.declare("tech_config", types=dict) + + def setup(self): + # Load in user config from input file + self.config = ShellTubeHXPerformanceModelConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), + additional_cls_name=self.__class__.__name__, + ) + + # Setup OpenMDAO inputs + self.add_input( + "Th_in_C", val=self.config.Th_in_C, units="C", desc="Hot fluid inlet temperature" + ) + self.add_input( + "Tc_in_C", val=self.config.Tc_in_C, units="C", desc="Cold fluid inlet temperature" + ) + self.add_input( + "m_dot_h_kg_s", + val=self.config.m_dot_h_kg_s, + units="kg/s", + desc="Hot fluid mass flow rate", + ) + self.add_input( + "m_dot_c_kg_s", + val=self.config.m_dot_c_kg_s, + units="kg/s", + desc="Cold fluid mass flow rate", + ) + + # Setup OpenMDAO outputs + self.add_output("Th_out_C", val=0.0, units="C", desc="Hot fluid outlet temperature") + self.add_output("Tc_out_C", val=0.0, units="C", desc="Cold fluid outlet temperature") + self.add_output("Q_total_kW", val=0.0, units="kW", desc="Total heat transfer rate") + self.add_output("epsilon", val=0.0, desc="Effectiveness of the heat exchanger") + self.add_output("NTU", val=0.0, desc="Number of transfer units") + self.add_output("C_r", val=0.0, desc="Heat capacity rate ratio") + self.add_output( + "U_global_W_m2K", val=0.0, units="W/m**2/K", desc="Overall heat transfer coefficient" + ) + self.add_output("dp_hot_Pa", val=0.0, units="Pa", desc="Pressure drop on the hot side") + self.add_output("dp_cold_Pa", val=0.0, units="Pa", desc="Pressure drop on the cold side") + self.add_output("pump_power_kW", val=0.0, units="kW", desc="Total pump power required") + self.add_output("S_gen_dot_W_per_K", val=0.0, units="W/K", desc="Entropy generation rate") + self.add_output("Ex_dest_dot_kW", val=0.0, units="kW", desc="Exergy destruction rate") + + def compute(self, inputs, outputs): + params = { + "Th_in": inputs["Th_in_C"][0], + "Tc_in": inputs["Tc_in_C"][0], + "m_dot_h": inputs["m_dot_h_kg_s"][0], + "m_dot_c": inputs["m_dot_c_kg_s"][0], + "N_tubes": self.config.N_tubes, + "N_passes": self.config.N_passes, + "L_tube": self.config.L_tube_m, + "D_o": self.config.D_o_m, + "t_wall": self.config.t_wall_m, + "D_shell": self.config.D_shell_m, + } + + res = hx_shell_tube_steady(params) + + outputs["Th_out_C"] = float(res["Th"][-1]) + outputs["Tc_out_C"] = float(res["Tc"][0]) + outputs["Q_total_kW"] = float(res["Q_total"]) / 1e3 + outputs["epsilon"] = float(res["epsilon"]) + outputs["NTU"] = float(res["NTU"]) + outputs["C_r"] = float(res["C_r"]) + outputs["U_global_W_m2K"] = float(res["U_global"]) + outputs["dp_hot_Pa"] = float(res["dp_tube_total"]) + outputs["dp_cold_Pa"] = float(res["dp_shell_total"]) + outputs["pump_power_kW"] = float(res["P_pump_total"]) / 1e3 + outputs["S_gen_dot_W_per_K"] = float(res["S_gen_dot"]) + outputs["Ex_dest_dot_kW"] = float(res["Ex_dest_dot"]) / 1e3 diff --git a/h2integrate/converters/heat/shell_tube_hx_cost_model.py b/h2integrate/converters/heat/shell_tube_hx_cost_model.py new file mode 100644 index 000000000..98ab609ad --- /dev/null +++ b/h2integrate/converters/heat/shell_tube_hx_cost_model.py @@ -0,0 +1,71 @@ +from attrs import field, define + +from h2integrate.core.utilities import merge_shared_inputs +from h2integrate.core.validators import gt_zero +from h2integrate.core.model_baseclasses import CostModelBaseClass, CostModelBaseConfig + + +@define(kw_only=True) +class ShellTubeHXCostModelConfig(CostModelBaseConfig): + """ + Configuration class for the ShellTubeHXCostModel. + + Args: + Q_ref (float|int): Reference heat transfer rate for cost scaling in W + C_ref (float|int): Reference capital cost for cost scaling in USD + exp_Q (float|int): Exponent for heat transfer cost scaling + """ + + Q_ref: float | int = field(validator=gt_zero) + C_ref: float | int = field(validator=gt_zero) + exp_Q: float | int = field(validator=gt_zero) + + +class ShellTubeHXCostModel(CostModelBaseClass): + """ + This is a very simple placeholder cost model: + + - Reference: 240,000 USD for a 1 MW HX with U ~ 1000 W/m²-K + - Scale CapEx ~ (Q / Q_ref)^0.8 + - OpEx = 4% of CapEx per year + """ + + def setup(self): + self.config = ShellTubeHXCostModelConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "cost"), + additional_cls_name=self.__class__.__name__, + ) + super().setup() + + self.add_input("Q_total_W", val=0.0, units="W", desc="Total heat transfer rate") + self.add_input( + "Q_ref", + val=self.config.Q_ref, + units="W", + desc="Reference heat transfer rate for cost scaling", + ) + self.add_input( + "C_ref", + val=self.config.C_ref, + units="USD", + desc="Reference capital cost for cost scaling", + ) + self.add_input( + "exp_Q", + val=self.config.exp_Q, + units="unitless", + desc="Exponent for heat transfer cost scaling", + ) + + def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): + Q_total_W = inputs["Q_total_W"][0] + Q_ref = inputs["Q_ref"][0] + C_ref = inputs["C_ref"][0] + exp_Q = inputs["exp_Q"][0] + + scale_Q = (Q_total_W / Q_ref) ** exp_Q if Q_total_W > 0 else 0.0 + capex = C_ref * scale_Q + opex = 0.04 * capex + + outputs["CapEx"] = capex + outputs["OpEx"] = opex diff --git a/h2integrate/converters/heat/test/test_etes.py b/h2integrate/converters/heat/test/test_etes.py new file mode 100644 index 000000000..3c3b5196e --- /dev/null +++ b/h2integrate/converters/heat/test/test_etes.py @@ -0,0 +1,212 @@ +import numpy as np +import pytest +import openmdao.api as om +from pytest import approx, fixture + +from h2integrate.converters.heat.etes import ETESPerformanceModel + + +@fixture +def petes_config(): + """P-ETES (decoupled) configuration using parameter values from Table 2.""" + performance_parameters = { + "etes_type": "P-ETES", + "S_TES_kWh": 1_000_000.0, # 1 GWh thermal storage + "S_ch_kW": 100_000.0, # 100 MW electric charging + "S_dis_kW": 50_000.0, # 50 MW thermal discharging + "eta_ch": 0.95, + "eta_dis": 0.73, + "f_loss": 0.0068, + "SOC_min": 0.0, + "SOC_init": 0.5, + "cost_year": 2026, + } + n = 24 + plant_config = { + "plant": {"simulation": {"n_timesteps": n, "dt": 3600}}, + } + tech_config = {"model_inputs": {"performance_parameters": performance_parameters}} + return tech_config, plant_config, n + + +@fixture +def retes_config(): + """R-ETES (integrated) configuration using parameter values from Table 3.""" + performance_parameters = { + "etes_type": "R-ETES", + "S_TES_kWh": 1_000_000.0, + "eta_ch": 0.98, + "eta_dis": 0.90, + "f_loss": 0.0068, + "f_ch_max": 0.3, + "f_dis_max": 0.2, + "SOC_min": 0.0, + "SOC_init": 0.5, + "cost_year": 2026, + } + n = 24 + plant_config = { + "plant": {"simulation": {"n_timesteps": n, "dt": 3600}}, + } + tech_config = {"model_inputs": {"performance_parameters": performance_parameters}} + return tech_config, plant_config, n + + +def _create_problem(tech_config, plant_config): + prob = om.Problem() + prob.model.add_subsystem( + "etes", + ETESPerformanceModel( + tech_config=tech_config, + plant_config=plant_config, + driver_config={}, + ), + promotes=["*"], + ) + prob.setup() + return prob + + +@pytest.mark.unit +class TestETESPerformanceModel: + def test_petes_discharge_meets_load(self, petes_config): + tech_config, plant_config, n = petes_config + prob = _create_problem(tech_config, plant_config) + + # No charging electricity, constant thermal load well below discharge cap + load = np.full(n, 10_000.0) # 10 MW_th + prob.set_val("electricity_in_kW", np.zeros(n)) + prob.set_val("heat_demand_kW", load) + prob.run_model() + + heat_out = prob.get_val("heat_out_kW", units="kW") + unmet = prob.get_val("unmet_heat_demand_kW", units="kW") + E_st = prob.get_val("E_st_kWh", units="kW*h") + + # Discharge should fully meet load since storage is half-full + assert heat_out == approx(load, rel=1e-6) + assert unmet == approx(np.zeros(n), abs=1e-6) + # Storage should monotonically decrease + assert np.all(np.diff(E_st) <= 0) + + def test_petes_charging_rate_limited(self, petes_config): + tech_config, plant_config, n = petes_config + prob = _create_problem(tech_config, plant_config) + + # Massive electricity input, no load. Charging should be capped at + # S_ch * eta_ch = 100_000 * 0.95 = 95_000 kW_th into storage. + prob.set_val("electricity_in_kW", np.full(n, 1e9)) + prob.set_val("heat_demand_kW", np.zeros(n)) + prob.run_model() + + Q_ch = prob.get_val("Q_ch_kW", units="kW") + S_ch_kW = tech_config["model_inputs"]["performance_parameters"]["S_ch_kW"] + eta_ch = tech_config["model_inputs"]["performance_parameters"]["eta_ch"] + cap = S_ch_kW * eta_ch + # Either capped by S_ch * eta_ch or by remaining headroom in storage + assert np.all(Q_ch <= cap + 1e-6) + # In early timesteps before storage fills, should hit the cap + assert Q_ch[0] == approx(cap, rel=1e-6) + + def test_petes_discharge_rate_limited(self, petes_config): + tech_config, plant_config, n = petes_config + # Start full so we have plenty of energy + tech_config["model_inputs"]["performance_parameters"]["SOC_init"] = 1.0 + prob = _create_problem(tech_config, plant_config) + + # Huge load, no charging + prob.set_val("electricity_in_kW", np.zeros(n)) + prob.set_val("heat_demand_kW", np.full(n, 1e9)) + prob.run_model() + + heat_out = prob.get_val("heat_out_kW", units="kW") + # Should be capped at S_dis_kW = 50_000 kW_th + S_dis_kW = tech_config["model_inputs"]["performance_parameters"]["S_dis_kW"] + assert heat_out[0] == approx(S_dis_kW, rel=1e-6) + + def test_petes_energy_balance(self, petes_config): + tech_config, plant_config, n = petes_config + prob = _create_problem(tech_config, plant_config) + + rng = np.random.default_rng(0) + elec_in = rng.uniform(0, 80_000, n) + load = rng.uniform(0, 30_000, n) + prob.set_val("electricity_in_kW", elec_in) + prob.set_val("heat_demand_kW", load) + prob.run_model() + + E_st = prob.get_val("E_st_kWh", units="kW*h") + Q_ch = prob.get_val("Q_ch_kW", units="kW") + Q_dis = prob.get_val("Q_dis_kW", units="kW") + Q_st_loss = prob.get_val("Q_st_loss_kW", units="kW") + dt = 1.0 # hour + + SOC_init = tech_config["model_inputs"]["performance_parameters"]["SOC_init"] + S_TES = tech_config["model_inputs"]["performance_parameters"]["S_TES_kWh"] + E_prev = SOC_init * S_TES + for t in range(n): + # E_st(t) = E_st(t-1) + (Q_ch - Q_dis - Q_st_loss) * dt + expected = E_prev + (Q_ch[t] - Q_dis[t] - Q_st_loss[t]) * dt + assert E_st[t] == approx(expected, rel=1e-9, abs=1e-6) + E_prev = E_st[t] + + def test_petes_loss_rates(self, petes_config): + tech_config, plant_config, n = petes_config + prob = _create_problem(tech_config, plant_config) + + prob.set_val("electricity_in_kW", np.full(n, 50_000.0)) + prob.set_val("heat_demand_kW", np.full(n, 20_000.0)) + prob.run_model() + + Q_ch = prob.get_val("Q_ch_kW", units="kW") + Q_dis = prob.get_val("Q_dis_kW", units="kW") + Q_ch_loss = prob.get_val("Q_ch_loss_kW", units="kW") + Q_dis_loss = prob.get_val("Q_dis_loss_kW", units="kW") + + eta_ch = tech_config["model_inputs"]["performance_parameters"]["eta_ch"] + eta_dis = tech_config["model_inputs"]["performance_parameters"]["eta_dis"] + + # Q_ch_loss = Q_ch * (1 - eta_ch) / eta_ch + assert Q_ch_loss == approx(Q_ch * (1 - eta_ch) / eta_ch, rel=1e-9) + # Q_dis_loss = Q_dis * (1 - eta_dis) + assert Q_dis_loss == approx(Q_dis * (1 - eta_dis), rel=1e-9) + + def test_retes_rate_caps(self, retes_config): + tech_config, plant_config, n = retes_config + # Start full to test discharge cap + tech_config["model_inputs"]["performance_parameters"]["SOC_init"] = 1.0 + prob = _create_problem(tech_config, plant_config) + + prob.set_val("electricity_in_kW", np.full(n, 1e9)) + prob.set_val("heat_demand_kW", np.full(n, 1e9)) + prob.run_model() + + Q_ch = prob.get_val("Q_ch_kW", units="kW") + Q_dis = prob.get_val("Q_dis_kW", units="kW") + params = tech_config["model_inputs"]["performance_parameters"] + S_TES = params["S_TES_kWh"] + # Rates capped at f * S_TES + assert np.all(Q_ch <= params["f_ch_max"] * S_TES + 1e-6) + assert np.all(Q_dis <= params["f_dis_max"] * S_TES + 1e-6) + + def test_soc_bounds(self, petes_config): + tech_config, plant_config, n = petes_config + tech_config["model_inputs"]["performance_parameters"]["SOC_min"] = 0.1 + tech_config["model_inputs"]["performance_parameters"]["SOC_init"] = 1.0 + prob = _create_problem(tech_config, plant_config) + + prob.set_val("electricity_in_kW", np.full(n, 1e9)) + prob.set_val("heat_demand_kW", np.full(n, 1e9)) + prob.run_model() + + E_st = prob.get_val("E_st_kWh", units="kW*h") + S_TES = tech_config["model_inputs"]["performance_parameters"]["S_TES_kWh"] + SOC_min = tech_config["model_inputs"]["performance_parameters"]["SOC_min"] + assert np.all(E_st >= SOC_min * S_TES - 1e-6) + assert np.all(E_st <= S_TES + 1e-6) + + def test_invalid_etes_type(self, petes_config): + tech_config, plant_config, _ = petes_config + tech_config["model_inputs"]["performance_parameters"]["etes_type"] = "BAD" + with pytest.raises(ValueError, match="etes_type"): + _create_problem(tech_config, plant_config) diff --git a/h2integrate/converters/heat/test/test_etes_cost_model.py b/h2integrate/converters/heat/test/test_etes_cost_model.py new file mode 100644 index 000000000..00d81eb61 --- /dev/null +++ b/h2integrate/converters/heat/test/test_etes_cost_model.py @@ -0,0 +1,114 @@ +import pytest +import openmdao.api as om +from pytest import approx, fixture + +from h2integrate.converters.heat.etes_cost_model import ETESCostModel + + +@fixture +def petes_cost_config(): + """P-ETES cost config with placeholder linear coefficients.""" + cost_parameters = { + "C_lin_TES": 5.0, # $/kWh_th + "C_const_TES": 1e6, # $ + "C_min_TES": 1.5, # $/kWh_th (from Table 2) + "C_lin_ch": 100.0, # $/kW_e + "C_const_ch": 5e5, + "C_lin_dis": 150.0, # $/kW_th + "C_const_dis": 3e5, + "opex_fraction": 0.04, + "cost_year": 2026, + } + tech_config = {"model_inputs": {"cost_parameters": cost_parameters}} + plant_config = {"plant": {"plant_life": 30, "simulation": {"dt": 3600}}} + return tech_config, plant_config + + +@fixture +def retes_cost_config(): + """R-ETES cost config: no separate charging/discharging unit costs.""" + cost_parameters = { + "C_lin_TES": 10.0, + "C_const_TES": 5e5, + "C_min_TES": 3.0, + "opex_fraction": 0.04, + "cost_year": 2026, + } + tech_config = {"model_inputs": {"cost_parameters": cost_parameters}} + plant_config = {"plant": {"plant_life": 30, "simulation": {"dt": 3600}}} + return tech_config, plant_config + + +def _create_problem(tech_config, plant_config): + prob = om.Problem() + prob.model.add_subsystem( + "etes_cost", + ETESCostModel( + tech_config=tech_config, + plant_config=plant_config, + driver_config={}, + ), + promotes=["*"], + ) + prob.setup() + return prob + + +@pytest.mark.unit +class TestETESCostModel: + def test_petes_cost_above_floor(self, petes_cost_config): + tech_config, plant_config = petes_cost_config + prob = _create_problem(tech_config, plant_config) + + S_TES = 1_000_000.0 # kWh + S_ch = 100_000.0 # kW + S_dis = 50_000.0 # kW + prob.set_val("S_TES_kWh", S_TES) + prob.set_val("S_ch_kW", S_ch) + prob.set_val("S_dis_kW", S_dis) + prob.run_model() + + params = tech_config["model_inputs"]["cost_parameters"] + # C_lin_TES * S_TES = 5e6, floor = 1.5 * 1e6 = 1.5e6 -> linear wins + c_tes = params["C_lin_TES"] * S_TES + params["C_const_TES"] + c_ch = params["C_lin_ch"] * S_ch + params["C_const_ch"] + c_dis = params["C_lin_dis"] * S_dis + params["C_const_dis"] + expected_capex = c_tes + c_ch + c_dis + expected_opex = params["opex_fraction"] * expected_capex + + assert prob.get_val("CapEx", units="USD")[0] == approx(expected_capex, rel=1e-9) + assert prob.get_val("OpEx", units="USD/year")[0] == approx(expected_opex, rel=1e-9) + + def test_petes_cost_floor_active(self, petes_cost_config): + tech_config, plant_config = petes_cost_config + # Force a very small linear coefficient so the floor is active + tech_config["model_inputs"]["cost_parameters"]["C_lin_TES"] = 0.1 + tech_config["model_inputs"]["cost_parameters"]["C_const_TES"] = 0.0 + prob = _create_problem(tech_config, plant_config) + + S_TES = 1_000_000.0 + prob.set_val("S_TES_kWh", S_TES) + prob.set_val("S_ch_kW", 0.0) + prob.set_val("S_dis_kW", 0.0) + prob.run_model() + + # C_min_TES * S_TES = 1.5 * 1e6 = 1.5e6 (vs linear 0.1 * 1e6 = 1e5) + params = tech_config["model_inputs"]["cost_parameters"] + expected = params["C_min_TES"] * S_TES + assert prob.get_val("CapEx", units="USD")[0] == approx(expected, rel=1e-9) + + def test_retes_no_ch_dis_costs(self, retes_cost_config): + tech_config, plant_config = retes_cost_config + prob = _create_problem(tech_config, plant_config) + + S_TES = 500_000.0 + prob.set_val("S_TES_kWh", S_TES) + # R-ETES leaves S_ch/S_dis at zero + prob.run_model() + + params = tech_config["model_inputs"]["cost_parameters"] + c_tes_linear = params["C_lin_TES"] * S_TES + params["C_const_TES"] + c_tes_floor = params["C_min_TES"] * S_TES + expected_capex = max(c_tes_linear, c_tes_floor) + + assert prob.get_val("CapEx", units="USD")[0] == approx(expected_capex, rel=1e-9) diff --git a/h2integrate/converters/heat/test/test_etes_milp.py b/h2integrate/converters/heat/test/test_etes_milp.py new file mode 100644 index 000000000..7f442adb7 --- /dev/null +++ b/h2integrate/converters/heat/test/test_etes_milp.py @@ -0,0 +1,257 @@ +import numpy as np +import pytest +from pytest import approx + +from h2integrate.converters.heat.etes_milp import ETESMILPConfig, solve_etes_milp + + +def _peak_price_profile(n=24): + """Cheap nights, expensive afternoons (peaks 12-20).""" + price = np.full(n, 0.02) + price[12:20] = 0.20 + return price + + +@pytest.mark.unit +class TestETESMILP: + def test_petes_sizing_basic(self): + n = 24 + load = np.full(n, 10_000.0) # 10 MW_th constant + price = _peak_price_profile(n) + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.0, + t_ch_min_h=2.0, + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + fixed_charge_rate=0.10, + opex_fraction=0.04, + cyclic=True, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + + # Sizes must be strictly positive and finite + assert res.S_TES_kWh > 0 + assert res.S_ch_kW > 0 + assert res.S_dis_kW > 0 + # Total demand over horizon = 10 MW * 24 h = 240 MWh + total_heat_delivered = float(np.sum(res.Q_dis_kW * cfg.eta_dis)) + assert total_heat_delivered == approx(np.sum(load), rel=1e-4) + # No unmet load + assert np.all(res.unmet_load_kW < 1e-6) + + def test_petes_uses_cheap_electricity(self): + n = 24 + load = np.full(n, 10_000.0) + price = _peak_price_profile(n) + # Fix sizes so the optimizer doesn't trade off S_ch cost vs. electricity + # cost. With ample charging capacity, charging should fully avoid the + # expensive hours. + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.0, + t_ch_min_h=2.0, + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + S_TES_fixed_kWh=500_000.0, + S_ch_fixed_kW=100_000.0, + S_dis_fixed_kW=20_000.0, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + + # Charging should fully avoid the expensive hours (price = 0.20) + cheap_hours = price < 0.05 + expensive_hours = ~cheap_hours + ch_in_cheap = float(np.sum(res.Q_ch_kW[cheap_hours])) + ch_in_expensive = float(np.sum(res.Q_ch_kW[expensive_hours])) + assert ch_in_expensive < 1e-3 + assert ch_in_cheap > 0 + + def test_storage_energy_balance(self): + n = 24 + load = np.full(n, 8_000.0) + price = _peak_price_profile(n) + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.01, + t_ch_min_h=2.0, + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + + dt = 1.0 + E_prev = cfg.SOC_init * res.S_TES_kWh + for t in range(n): + expected = E_prev * (1.0 - cfg.f_loss) + (res.Q_ch_kW[t] - res.Q_dis_kW[t]) * dt + assert res.E_st_kWh[t] == approx(expected, rel=1e-5, abs=1e-3) + E_prev = res.E_st_kWh[t] + + def test_cyclic_constraint(self): + n = 12 + load = np.full(n, 5_000.0) + price = np.tile([0.02, 0.20], n // 2) + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.0, + t_ch_min_h=2.0, + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + cyclic=True, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + # Final E_st == initial E_st (= SOC_init * S_TES) + E0 = cfg.SOC_init * res.S_TES_kWh + assert res.E_st_kWh[-1] == approx(E0, rel=1e-5, abs=1e-3) + + def test_retes_rate_coupling(self): + n = 24 + load = np.full(n, 10_000.0) + price = _peak_price_profile(n) + cfg = ETESMILPConfig( + etes_type="R-ETES", + eta_ch=0.98, + eta_dis=0.90, + f_loss=0.0068, + f_ch_max=0.3, + f_dis_max=0.2, + C_lin_TES=10.0, + C_min_TES=3.0, + fixed_charge_rate=0.10, + opex_fraction=0.04, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + + # In R-ETES, S_ch and S_dis are not sized (= 0); charging/discharging + # rates are bounded by f * S_TES. + assert res.S_ch_kW == approx(0.0, abs=1e-6) + assert res.S_dis_kW == approx(0.0, abs=1e-6) + assert np.all(res.Q_ch_kW <= cfg.f_ch_max * res.S_TES_kWh + 1e-4) + assert np.all(res.Q_dis_kW <= cfg.f_dis_max * res.S_TES_kWh + 1e-4) + + def test_minimum_charging_time(self): + n = 24 + load = np.full(n, 10_000.0) + price = _peak_price_profile(n) + t_min = 6.0 + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.0, + t_ch_min_h=t_min, + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + # S_ch * eta_ch * t_ch_min <= S_TES + assert res.S_ch_kW * cfg.eta_ch * t_min <= res.S_TES_kWh + 1e-3 + + def test_dispatch_only_mode(self): + """Fix sizes and only optimize dispatch.""" + n = 24 + load = np.full(n, 8_000.0) + price = _peak_price_profile(n) + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.0, + t_ch_min_h=2.0, + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + S_TES_fixed_kWh=200_000.0, + S_ch_fixed_kW=50_000.0, + S_dis_fixed_kW=20_000.0, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + assert res.S_TES_kWh == approx(200_000.0, rel=1e-6) + assert res.S_ch_kW == approx(50_000.0, rel=1e-6) + assert res.S_dis_kW == approx(20_000.0, rel=1e-6) + # Charging unit not exceeded + assert np.all(res.Q_ch_kW <= res.S_ch_kW * cfg.eta_ch + 1e-4) + # Discharging unit not exceeded + assert np.all(res.Q_dis_kW * cfg.eta_dis <= res.S_dis_kW + 1e-4) + + def test_unmet_load_allowed(self): + """If sizes are too small and unmet allowed, model should drop load + when penalty < grid cost.""" + n = 24 + load = np.full(n, 100_000.0) # very large load + price = np.full(n, 0.50) # expensive grid + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.0, + t_ch_min_h=0.0, # disable to avoid sizing coupling + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + S_TES_fixed_kWh=1.0, # essentially no storage + S_ch_fixed_kW=1.0, + S_dis_fixed_kW=1.0, + allow_unmet_load=True, + unmet_load_penalty=0.10, # cheaper to shed than charge from grid + cyclic=False, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + # Should have substantial unmet load + assert np.sum(res.unmet_load_kW) > 0 + + def test_cost_breakdown_consistency(self): + n = 24 + load = np.full(n, 5_000.0) + price = _peak_price_profile(n) + cfg = ETESMILPConfig( + etes_type="P-ETES", + eta_ch=0.95, + eta_dis=0.73, + f_loss=0.0, + t_ch_min_h=2.0, + C_lin_TES=5.0, + C_min_TES=1.5, + C_lin_ch=100.0, + C_lin_dis=150.0, + fixed_charge_rate=0.10, + opex_fraction=0.04, + ) + res = solve_etes_milp(cfg, price, load, dt_h=1.0) + + # Total annualized cost = annualized capex + opex + electricity + assert res.total_annualized_cost_USD == approx( + res.annualized_capex_USD + res.annual_opex_USD + res.annual_electricity_cost_USD, + rel=1e-6, + ) + # Objective value matches total (no unmet load => no penalty) + assert res.objective_value == approx(res.total_annualized_cost_USD, rel=1e-6) + + def test_invalid_etes_type_raises(self): + with pytest.raises(ValueError, match="etes_type"): + ETESMILPConfig( + etes_type="BAD", + eta_ch=0.9, + eta_dis=0.7, + f_loss=0.0, + ) diff --git a/h2integrate/converters/heat/test/test_shell_tube_hx.py b/h2integrate/converters/heat/test/test_shell_tube_hx.py new file mode 100644 index 000000000..eef1e45ae --- /dev/null +++ b/h2integrate/converters/heat/test/test_shell_tube_hx.py @@ -0,0 +1,90 @@ +import pytest +import openmdao.api as om +from pytest import approx, fixture + +from h2integrate.converters.heat.shell_tube_hx import ShellTubeHXPerformanceModel + + +@fixture +def shell_tube_hx_config(): + performance_parameterss = { + "Th_in_C": 90.0, + "Tc_in_C": 30.0, + "m_dot_h_kg_s": 18.0, + "m_dot_c_kg_s": 40.0, + "N_tubes": 192, + "N_passes": 2, + "L_tube_m": 4.9, + "D_o_m": 0.01905, + "t_wall_m": 0.002, + "D_shell_m": 0.591, + "cost_year": 2024, + } + + tech_config = { + "model_inputs": { + "performance_parameters": performance_parameterss, + } + } + return tech_config + + +@pytest.mark.unit +class TestShellTubeHXPerformanceModel: + def _create_problem(self, config): + prob = om.Problem() + prob.model.add_subsystem( + "shell_tube_hx", + ShellTubeHXPerformanceModel(tech_config=config), + promotes=["*"], + ) + prob.setup() + return prob + + def test_hx_performance_calculation(self, shell_tube_hx_config): + prob = self._create_problem(shell_tube_hx_config) + + prob.run_model() + + C_r = prob.get_val("shell_tube_hx.C_r") + Ex_dest_dot_kW = prob.get_val("shell_tube_hx.Ex_dest_dot_kW", units="kW") + NTU = prob.get_val("shell_tube_hx.NTU") + Q_total_kW = prob.get_val("shell_tube_hx.Q_total_kW", units="kW") + S_gen_dot_W_per_K = prob.get_val("shell_tube_hx.S_gen_dot_W_per_K", units="W/K") + Th_out_C = prob.get_val("shell_tube_hx.Th_out_C", units="C") + Tc_out_C = prob.get_val("shell_tube_hx.Tc_out_C", units="C") + U_global_W_m2K = prob.get_val("shell_tube_hx.U_global_W_m2K", units="W/m**2/K") + dp_cold_Pa = prob.get_val("shell_tube_hx.dp_cold_Pa", units="Pa") + dp_hot_Pa = prob.get_val("shell_tube_hx.dp_hot_Pa", units="Pa") + epsilon = prob.get_val("shell_tube_hx.epsilon") + pump_power_kW = prob.get_val("shell_tube_hx.pump_power_kW", units="kW") + + # Expected values can be adjusted based on known results or calculations + expected_C_r = 0.4527329816950921 + expected_Ex_dest_dot_kW = 248.78629912498585 + expected_NTU = 0.9142424832326593 + expected_Q_total_kW = 2462.235075109509 + expected_S_gen_dot_W_per_K = 834.433335988549 + expected_Tc_out_C = 44.72803509154916 + expected_Th_out_C = 57.381403669025964 + expected_U_global_W_m2K = 1229.0775390792805 + expected_dp_cold_Pa = 201.87162586061777 + expected_dp_hot_Pa = 10251.659491792187 + expected_epsilon = 0.5421484468743297 + expected_pump_power_kW = 0.28157427036731625 + + rel = 1e-5 + assert Q_total_kW == approx(expected_Q_total_kW, rel=rel) + assert U_global_W_m2K == approx(expected_U_global_W_m2K, rel=rel) + assert C_r == approx(expected_C_r, rel=rel) + assert Ex_dest_dot_kW == approx(expected_Ex_dest_dot_kW, rel=rel) + assert NTU == approx(expected_NTU, rel=rel) + assert Q_total_kW == approx(expected_Q_total_kW, rel=rel) + assert S_gen_dot_W_per_K == approx(expected_S_gen_dot_W_per_K, rel=rel) + assert Tc_out_C == approx(expected_Tc_out_C, rel=rel) + assert Th_out_C == approx(expected_Th_out_C, rel=rel) + assert U_global_W_m2K == approx(expected_U_global_W_m2K, rel=rel) + assert dp_cold_Pa == approx(expected_dp_cold_Pa, rel=rel) + assert dp_hot_Pa == approx(expected_dp_hot_Pa, rel=rel) + assert epsilon == approx(expected_epsilon, rel=rel) + assert pump_power_kW == approx(expected_pump_power_kW, rel=rel) diff --git a/h2integrate/converters/heat/test/test_shell_tube_hx_cost_model.py b/h2integrate/converters/heat/test/test_shell_tube_hx_cost_model.py new file mode 100644 index 000000000..f30d5f708 --- /dev/null +++ b/h2integrate/converters/heat/test/test_shell_tube_hx_cost_model.py @@ -0,0 +1,67 @@ +import pytest +import openmdao.api as om +from pytest import approx, fixture + +from h2integrate.converters.heat.shell_tube_hx_cost_model import ShellTubeHXCostModel + + +@fixture +def shell_tube_hx_cost_model_config(): + cost_parameters = { + "Q_ref": 1.0e6, + "C_ref": 2.4e5, + "exp_Q": 0.8, + "cost_year": 2024, + } + + tech_config = { + "model_inputs": { + "cost_parameters": cost_parameters, + } + } + + plant_config = { + "plant": { + "plant_life": 30, + "simulation": { + "n_timesteps": 8760, + "dt": 3600, + }, + } + } + return tech_config, plant_config + + +@pytest.mark.unit +class TestShellTubeHXCostModel: + def _create_problem(self, config): + prob = om.Problem() + prob.model.add_subsystem( + "shell_tube_hx_cost", + ShellTubeHXCostModel( + tech_config=config[0], + plant_config=config[1], + driver_config={}, + ), + promotes=["*"], + ) + prob.setup() + return prob + + def test_hx_cost_calculation(self, shell_tube_hx_cost_model_config): + prob = self._create_problem(shell_tube_hx_cost_model_config) + + # Set input value for total heat transfer rate + prob.set_val("Q_total_W", 5.0e6) # 5 MW + + prob.run_model() + + CapEx_USD = prob.get_val("shell_tube_hx_cost.CapEx", units="USD") + OpEx_USD_per_year = prob.get_val("shell_tube_hx_cost.OpEx", units="USD/year") + + expected_CapEx_USD = 240000.0 * (5.0e6 / 1.0e6) ** 0.8 + expected_OpEx_USD_per_year = 0.04 * expected_CapEx_USD + + rel = 1e-5 + assert CapEx_USD == approx(expected_CapEx_USD, rel=rel) + assert OpEx_USD_per_year == approx(expected_OpEx_USD_per_year, rel=rel)