diff --git a/config/config.default.yaml b/config/config.default.yaml index 0af34734..c9ee439b 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -418,6 +418,15 @@ sector: 2045: 0.8 2050: 1.0 district_heating_loss: 0.15 + # check these numbers! + forward_temperature: 90 #C + return_temperature: 60 #C + heat_source_cooling: 6 #K + heat_pump_cop_approximation: + refrigerant: ammonia + heat_exchanger_pinch_point_temperature_difference: 5 #K + isentropic_compressor_efficiency: 0.8 + heat_loss: 0.0 cluster_heat_buses: true heat_demand_cutout: default bev_dsm_restriction_value: 0.75 @@ -500,7 +509,7 @@ sector: aviation_demand_factor: 1. HVC_demand_factor: 1. time_dep_hp_cop: true - heat_pump_sink_T: 55. + heat_pump_sink_T_individual_heating: 55. reduce_space_heat_exogenously: true reduce_space_heat_exogenously_factor: 2020: 0.10 # this results in a space heat demand reduction of 10% diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 9986bceb..41c857b3 100644 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -217,13 +217,19 @@ rule build_temperature_profiles: rule build_cop_profiles: params: - heat_pump_sink_T=config_provider("sector", "heat_pump_sink_T"), + heat_pump_sink_T_individual_heating=config_provider("sector", "heat_pump_sink_T_individual_heating"), + forward_temperature_district_heating=config_provider("sector", "district_heating", "forward_temperature"), + return_temperature_district_heating=config_provider("sector", "district_heating", "return_temperature"), + heat_source_cooling_district_heating=config_provider("sector", "district_heating", "heat_source_cooling"), + heat_pump_cop_approximation=config_provider("sector", "district_heating", "heat_pump_cop_approximation"), input: temp_soil_total=resources("temp_soil_total_elec_s{simpl}_{clusters}.nc"), temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"), output: cop_soil_individual_heating=resources("cop_soil_individual_heating_elec_s{simpl}_{clusters}.nc"), cop_air_individual_heating=resources("cop_air_individual_heating_elec_s{simpl}_{clusters}.nc"), + cop_air_district_heating=resources("cop_air_district_heating_elec_s{simpl}_{clusters}.nc"), + cop_soil_district_heating=resources("cop_soil_district_heating_elec_s{simpl}_{clusters}.nc"), resources: mem_mb=20000, log: @@ -1023,6 +1029,8 @@ rule prepare_sector_network: temp_air_urban=resources("temp_air_urban_elec_s{simpl}_{clusters}.nc"), cop_soil_individual_heating=resources("cop_soil_individual_heating_elec_s{simpl}_{clusters}.nc"), cop_air_individual_heating=resources("cop_air_individual_heating_elec_s{simpl}_{clusters}.nc"), + cop_air_district_heating=resources("cop_air_district_heating_elec_s{simpl}_{clusters}.nc"), + cop_soil_district_heating=resources("cop_soil_district_heating_elec_s{simpl}_{clusters}.nc"), solar_thermal_total=lambda w: ( resources("solar_thermal_total_elec_s{simpl}_{clusters}.nc") if config_provider("sector", "solar_thermal")(w) diff --git a/scripts/build_cop_profiles.py b/scripts/build_cop_profiles.py index 104737b8..6d7050ba 100644 --- a/scripts/build_cop_profiles.py +++ b/scripts/build_cop_profiles.py @@ -42,11 +42,14 @@ References [1] Staffell et al., Energy & Environmental Science 11 (2012): A review of domestic heat pumps, https://doi.org/10.1039/C2EE22653G. """ +from typing import Union +from enum import Enum import xarray as xr +import numpy as np from _helpers import set_scenario_config -def coefficient_of_performance(delta_T, source="air"): +def coefficient_of_performance_individual_heating(delta_T, source="air"): if source == "air": return 6.81 - 0.121 * delta_T + 0.000630 * delta_T**2 elif source == "soil": @@ -55,6 +58,301 @@ def coefficient_of_performance(delta_T, source="air"): raise NotImplementedError("'source' must be one of ['air', 'soil']") +def celsius_to_kelvin(t_celsius: Union[float, xr.DataArray, np.array]) -> Union[float, xr.DataArray, np.array]: + if (np.asarray(t_celsius) > 200).any(): + raise ValueError("t_celsius > 200. Are you sure you are using the right units?") + return t_celsius + 273.15 + + +def logarithmic_mean(t_hot: Union[float, xr.DataArray, np.ndarray], t_cold: Union[float, xr.DataArray, np.ndarray]) -> Union[float, xr.DataArray, np.ndarray]: + if (np.asarray(t_hot <= t_cold)).any(): + raise ValueError("t_hot must be greater than t_cold") + return (t_hot - t_cold) / np.log(t_hot / t_cold) + +class CopDistrictHeating: + + def __init__( + self, + forward_temperature_celsius: Union[xr.DataArray, np.array], + source_inlet_temperature_celsius: Union[xr.DataArray, np.array], + return_temperature_celsius: Union[xr.DataArray, np.array], + source_outlet_temperature_celsius: Union[xr.DataArray, np.array], + delta_t_pinch_point: float = 5, + isentropic_compressor_efficiency: float = 0.8, + heat_loss: float = 0.0, + ) -> None: + """ + Initialize the COPProfileBuilder object. + + Parameters: + ---------- + forward_temperature_celsius : Union[xr.DataArray, np.array] + The forward temperature in Celsius. + return_temperature_celsius : Union[xr.DataArray, np.array] + The return temperature in Celsius. + source_inlet_temperature_celsius : Union[xr.DataArray, np.array] + The source inlet temperature in Celsius. + source_outlet_temperature_celsius : Union[xr.DataArray, np.array] + The source outlet temperature in Celsius. + delta_t_pinch_point : float, optional + The pinch point temperature difference, by default 5. + isentropic_compressor_efficiency : float, optional + The isentropic compressor efficiency, by default 0.8. + heat_loss : float, optional + The heat loss, by default 0.0. + """ + self.t_source_in = celsius_to_kelvin(source_inlet_temperature_celsius) + self.t_sink_out = celsius_to_kelvin(forward_temperature_celsius) + + self.t_sink_in = celsius_to_kelvin(return_temperature_celsius) + self.t_source_out = celsius_to_kelvin(source_outlet_temperature_celsius) + + self.isentropic_efficiency_compressor = isentropic_compressor_efficiency + self.heat_loss = heat_loss + self.delta_t_pinch = delta_t_pinch_point + + def cop(self) -> Union[xr.DataArray, np.array]: + """ + Calculate the coefficient of performance (COP) for the system. + + Returns: + Union[xr.DataArray, np.array]: The calculated COP values. + """ + return ( + self.ideal_lorenz_cop + * ( + ( + 1 + + (self.delta_t_refrigerant_sink + self.delta_t_pinch) + / self.t_sink_mean + ) + / ( + 1 + + ( + self.delta_t_refrigerant_sink + + self.delta_t_refrigerant_source + + 2 * self.delta_t_pinch + ) + / self.delta_t_lift + ) + ) + * self.isentropic_efficiency_compressor + * (1 - self.ratio_evaporation_compression_work) + + 1 + - self.isentropic_efficiency_compressor + - self.heat_loss + ) + + @property + def t_sink_mean(self) -> Union[xr.DataArray, np.array]: + """ + Calculate the logarithmic mean temperature difference between the cold and hot sinks. + + Returns + ------- + Union[xr.DataArray, np.array] + The mean temperature difference. + """ + return logarithmic_mean(t_cold=self.t_sink_in, t_hot=self.t_sink_out) + + @property + def t_source_mean(self) -> Union[xr.DataArray, np.array]: + """ + Calculate the logarithmic mean temperature of the heat source. + + Returns + ------- + Union[xr.DataArray, np.array] + The mean temperature of the heat source. + """ + return logarithmic_mean(t_hot=self.t_source_in, t_cold=self.t_source_out) + + @property + def delta_t_lift(self) -> Union[xr.DataArray, np.array]: + """ + Calculate the temperature lift as the difference between the logarithmic sink and source temperatures. + + Returns + ------- + Union[xr.DataArray, np.array] + The temperature difference between the sink and source. + """ + return self.t_sink_mean - self.t_source_mean + + @property + def ideal_lorenz_cop(self) -> Union[xr.DataArray, np.array]: + """ + Ideal Lorenz coefficient of performance (COP). + + The ideal Lorenz COP is calculated as the ratio of the mean sink temperature + to the lift temperature difference. + + Returns + ------- + np.array + The ideal Lorenz COP. + + """ + return self.t_sink_mean / self.delta_t_lift + + @property + def delta_t_refrigerant_source(self) -> Union[xr.DataArray, np.array]: + """ + Calculate the temperature difference between the refrigerant source inlet and outlet. + + Returns + ------- + Union[xr.DataArray, np.array] + The temperature difference between the refrigerant source inlet and outlet. + """ + return self._approximate_delta_t_refrigerant_source( + delta_t_source=self.t_source_in - self.t_source_out + ) + + @property + def delta_t_refrigerant_sink(self) -> Union[xr.DataArray, np.array]: + """ + Temperature difference between the refrigerant and the sink based on approximation. + + Returns + ------- + Union[xr.DataArray, np.array] + The temperature difference between the refrigerant and the sink. + """ + return self._approximate_delta_t_refrigerant_sink() + + @property + def ratio_evaporation_compression_work(self) -> Union[xr.DataArray, np.array]: + """ + Calculate the ratio of evaporation to compression work based on approximation. + + Returns + ------- + Union[xr.DataArray, np.array] + The calculated ratio of evaporation to compression work. + """ + return self._ratio_evaporation_compression_work_approximation() + + @property + def delta_t_sink(self) -> Union[xr.DataArray, np.array]: + """ + Calculate the temperature difference at the sink. + + Returns + ------- + Union[xr.DataArray, np.array] + The temperature difference at the sink. + """ + return self.t_sink_out - self.t_sink_in + + def _approximate_delta_t_refrigerant_source( + self, delta_t_source: Union[xr.DataArray, np.array] + ) -> Union[xr.DataArray, np.array]: + """ + Approximates the temperature difference between the refrigerant and the source. + + Parameters + ---------- + delta_t_source : Union[xr.DataArray, np.array] + The temperature difference for the refrigerant source. + + Returns + ------- + Union[xr.DataArray, np.array] + The approximate temperature difference for the refrigerant source. + """ + return delta_t_source / 2 + + def _approximate_delta_t_refrigerant_sink( + self, + refrigerant: str = "ammonia", + a: float = {"ammonia": 0.2, "isobutane": -0.0011}, + b: float = {"ammonia": 0.2, "isobutane": 0.3}, + c: float = {"ammonia": 0.016, "isobutane": 2.4}, + ) -> Union[xr.DataArray, np.array]: + """ + Approximates the temperature difference at the refrigerant sink. + + Parameters: + ---------- + refrigerant : str, optional + The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'. + a : float, optional + Coefficient for the temperature difference between the sink and source, default is 0.2. + b : float, optional + Coefficient for the temperature difference at the sink, default is 0.2. + c : float, optional + Constant term, default is 0.016. + + Returns: + ------- + Union[xr.DataArray, np.array] + The approximate temperature difference at the refrigerant sink. + + Notes: + ------ + This function assumes ammonia as the refrigerant. + + The approximate temperature difference at the refrigerant sink is calculated using the following formula: + a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c + + """ + if refrigerant not in a.keys(): + raise ValueError( + f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}" + ) + return ( + a[refrigerant] + * (self.t_sink_out - self.t_source_out + 2 * self.delta_t_pinch) + + b[refrigerant] * self.delta_t_sink + + c[refrigerant] + ) + + def _ratio_evaporation_compression_work_approximation( + self, + refrigerant: str = "ammonia", + a: float = {"ammonia": 0.0014, "isobutane": 0.0035}, + b: float = {"ammonia": -0.0015, "isobutane": -0.0033}, + c: float = {"ammonia": 0.039, "isobutane": 0.053}, + ) -> Union[xr.DataArray, np.array]: + """ + Calculate the ratio of evaporation to compression work approximation. + + Parameters: + ---------- + refrigerant : str, optional + The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'. + a : float, optional + Coefficient 'a' in the approximation equation. Default is 0.0014. + b : float, optional + Coefficient 'b' in the approximation equation. Default is -0.0015. + c : float, optional + Coefficient 'c' in the approximation equation. Default is 0.039. + + Returns: + ------- + Union[xr.DataArray, np.array] + The calculated ratio of evaporation to compression work. + + Notes: + ------ + This function assumes ammonia as the refrigerant. + + The approximation equation used is: + ratio = a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c + """ + if refrigerant not in a.keys(): + raise ValueError( + f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}" + ) + return ( + a[refrigerant] + * (self.t_sink_out - self.t_source_out + 2 * self.delta_t_pinch) + + b[refrigerant] * self.delta_t_sink + + c[refrigerant] + ) + + if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -70,8 +368,18 @@ if __name__ == "__main__": for source in ["air", "soil"]: source_T = xr.open_dataarray(snakemake.input[f"temp_{source}_total"]) - delta_T = snakemake.params.heat_pump_sink_T - source_T + delta_T = snakemake.params.heat_pump_sink_T_individual_heating - source_T - cop = coefficient_of_performance(delta_T, source) + cop_individual_heating = coefficient_of_performance_individual_heating(delta_T, source) + cop_individual_heating.to_netcdf(snakemake.output[f"cop_{source}_individual_heating"]) - cop.to_netcdf(snakemake.output[f"cop_{source}_individual_heating"]) + cop_district_heating = CopDistrictHeating( + forward_temperature_celsius=snakemake.params.forward_temperature_district_heating, + return_temperature_celsius=snakemake.params.return_temperature_district_heating, + source_inlet_temperature_celsius=source_T, + source_outlet_temperature_celsius=source_T - snakemake.params.heat_source_cooling_district_heating, + ).cop() + + cop_district_heating.to_netcdf(snakemake.output[f"cop_{source}_district_heating"]) + +