From 61c843bb938e6c5a27eeca447861e19bd9c6ca68 Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 13:42:05 +0100 Subject: [PATCH 01/36] add get name --- vclibpy/datamodels.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/vclibpy/datamodels.py b/vclibpy/datamodels.py index 969ac4c..c2da3fe 100644 --- a/vclibpy/datamodels.py +++ b/vclibpy/datamodels.py @@ -132,6 +132,10 @@ def convert_to_str_value_format(self, with_unit_and_description: bool) -> Dict[s return {f"{k} in {v.unit} ({v.description})": v.value for k, v in self.items()} return {k: v.value for k, v in self.items()} + def get_name(self): + return ";".join([k + "=" + str(round(v.value, 3)).replace(".", "_") + for k, v in self.get_variables().items()]) + class FlowsheetState(VariableContainer): """ From 321c99e34628c3a059b5aba599a3e06eeb77d37d Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 13:42:30 +0100 Subject: [PATCH 02/36] fix ntu issue for zero --- vclibpy/components/heat_exchangers/moving_boundary_ntu.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py index 19d29cd..d0ee85b 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py @@ -83,8 +83,11 @@ def iterate_area(self, dT_max, alpha_pri, alpha_sec, Q) -> float: area = 0.0 while True: - NTU = self.calc_NTU(area, k, m_flow_cp_min) - eps = self.calc_eps(R, NTU) + if self.flow_type == "cross" and area == 0.0: + eps = 0.0 + else: + NTU = self.calc_NTU(area, k, m_flow_cp_min) + eps = self.calc_eps(R, NTU) if eps >= eps_necessary: if _step <= _accuracy: break From 0fb86758cc1a6eb7fcfcc6c4eeeccba17f2741d5 Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 13:43:14 +0100 Subject: [PATCH 03/36] add VDI plate heat exchanger --- .../heat_transfer/vdi_atlas_plate_he.py | 126 ++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 vclibpy/components/heat_exchangers/heat_transfer/vdi_atlas_plate_he.py diff --git a/vclibpy/components/heat_exchangers/heat_transfer/vdi_atlas_plate_he.py b/vclibpy/components/heat_exchangers/heat_transfer/vdi_atlas_plate_he.py new file mode 100644 index 0000000..da28398 --- /dev/null +++ b/vclibpy/components/heat_exchangers/heat_transfer/vdi_atlas_plate_he.py @@ -0,0 +1,126 @@ +import logging +from dataclasses import dataclass + +import numpy as np + +from .heat_transfer import HeatTransfer +from vclibpy.media import TransportProperties + +logger = logging.getLogger(__name__) + + +@dataclass +class PlateHeatExchangerGeometry: + """ + Geometry of a fin and tube heat exchanger with two rows of pipes in a shifted arrangement. + + Source: VDI-Wärmeatlas, Druckverlust und Wärmeübergang in Plattenwärmeübertragern, 11. Auflage, S.1687 + + + """ + wall_thickness: float # thickness of a plate + lambda_w: float # thermal conductivity of plates + amplitude: float # amplitude of a wavy plate + wave_length: float # wavelength of a wavy plate + phi: float # embossing angle + width: float # width of a plate + height: float # height of a plate + n_plates: float # number of plates + + @property + def X(self) -> float: + """return the wave number""" + return 2 * np.pi * self.amplitude / self.wave_length + + @property + def enlargement_factor(self) -> float: + """return surface enlargement factor""" + return (1 + np.sqrt(1 + self.X ** 2) + 4 * np.sqrt(1 + self.X ** 2 / 2)) / 6 + + @property + def d_h(self) -> float: + """return hydraulic diameter""" + return 4 * self.amplitude / self.enlargement_factor + + @property + def A(self) -> float: + """return plate surface""" + return self.height * self.width * self.enlargement_factor * self.n_plates + + +class VDIAtlasPlateHeatTransfer(HeatTransfer): + """ + VDI-Atlas based heat transfer estimation. + Check `AirToWallTransfer` for further argument options. + + This class assumes two pipes with a shifted arrangement. + + Args: + A_cross (float): Free cross-sectional area. + characteristic_length (float): Characteristic length d_a. + geometry_parameters (AirSourceHeatExchangerGeometry): + Geometry parameters for heat exchanger according to VDI-Atlas + """ + + def __init__(self, + geometry_parameters: PlateHeatExchangerGeometry): + super().__init__() + self.geometry_parameters = geometry_parameters + + def calc(self, transport_properties: TransportProperties, m_flow: float) -> float: + """ + Calculate heat transfer coefficient from refrigerant to the wall of the heat exchanger. + + The flow is assumed to be always turbulent and is based on a calibrated + Nusselt correlation. + + Args: + transport_properties (TransportProperties): Transport properties of the fluid. + m_flow (float): Mass flow rate of the fluid. + + Returns: + float: Heat transfer coefficient from refrigerant to HE in W/(m^2*K). + """ + Re = self.calc_reynolds( + dynamic_viscosity=transport_properties.dyn_vis, + m_flow=m_flow + ) + Nu = self.calc_turbulent_plate_nusselt(Re, transport_properties.Pr) + return Nu * transport_properties.lam / self.geometry_parameters.d_h + + def calc_reynolds(self, dynamic_viscosity: float, m_flow: float) -> float: + """ + Calculate Reynolds number for flow between two plates. + + Args: + dynamic_viscosity (float): Dynamic viscosity of the fluid. + m_flow (float): Mass flow rate. + width (float): Characteristic length (e.g., width) of the plate. + enlargement_factor (float): surface enlargement factor. + + Returns: + float: Reynolds number. + + """ + + return 2 * m_flow / (self.geometry_parameters.enlargement_factor * self.geometry_parameters.width * dynamic_viscosity) + + def calc_turbulent_plate_nusselt(self, Re: float, Pr: float) -> float: + """ + Calculate Nusselt Number based on Nusselt correlation VDI Atlas. + + Args: + Re (float): Reynolds number. + Pr (float): Prandtl number. + Zeta (float): pressure drop coefficient. + eta_by_eta_w (float): dynamic viscosity temperature correction factor. + + Returns: + float: Apparent heat transfer coefficient. + """ + c_q = 0.122 + q = 0.374 + eta_by_eta_w = 1 + Zeta = 1 + + return c_q * Pr ** (1/3) * eta_by_eta_w ** (1/6) * (Zeta * Re ** 2 * np.sin(2 * self.geometry_parameters.phi * np.pi / 180)) ** q From 26dd5db6c76c35791dc09d50bb7053aa2d46ab24 Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 13:44:12 +0100 Subject: [PATCH 04/36] Add first version for moving boundary based on logTm method --- .../heat_exchangers/heat_exchanger.py | 24 +- .../heat_exchangers/moving_boundary_Tm.py | 316 ++++++++++++++++++ vclibpy/components/heat_exchangers/ntu.py | 23 +- 3 files changed, 340 insertions(+), 23 deletions(-) create mode 100644 vclibpy/components/heat_exchangers/moving_boundary_Tm.py diff --git a/vclibpy/components/heat_exchangers/heat_exchanger.py b/vclibpy/components/heat_exchangers/heat_exchanger.py index 610ec8d..579ed49 100644 --- a/vclibpy/components/heat_exchangers/heat_exchanger.py +++ b/vclibpy/components/heat_exchangers/heat_exchanger.py @@ -37,11 +37,13 @@ def __init__( gas_heat_transfer: HeatTransfer, liquid_heat_transfer: HeatTransfer, two_phase_heat_transfer: TwoPhaseHeatTransfer, - secondary_medium: str + secondary_medium: str, + ratio_outer_to_inner_area: float = 1, ): super().__init__() self.A = A self.secondary_medium = secondary_medium.lower() + self.ratio_outer_to_inner_area = ratio_outer_to_inner_area self._wall_heat_transfer = wall_heat_transfer self._secondary_heat_transfer = secondary_heat_transfer @@ -164,6 +166,26 @@ def calc_alpha_secondary(self, transport_properties) -> float: m_flow=self.m_flow_secondary ) + def calc_k(self, alpha_pri: float, alpha_sec: float) -> float: + """ + Calculate the overall heat transfer coefficient (k) of the heat exchanger. + + Args: + alpha_pri (float): Heat transfer coefficient for the primary medium. + alpha_sec (float): Heat transfer coefficient for the secondary medium. + + Returns: + float: Overall heat transfer coefficient (k). + """ + k_wall = self.calc_wall_heat_transfer() + k = (1 / ( + (1 / alpha_pri) * self.ratio_outer_to_inner_area + + (1 / k_wall) * self.ratio_outer_to_inner_area + + (1 / alpha_sec) + ) + ) + return k + def calc_wall_heat_transfer(self) -> float: """ Calculate the heat transfer coefficient inside the wall. diff --git a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py new file mode 100644 index 0000000..1ebec10 --- /dev/null +++ b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py @@ -0,0 +1,316 @@ +import logging + +import numpy as np +from vclibpy.datamodels import FlowsheetState, Inputs +from vclibpy.components.heat_exchangers import HeatExchanger +from vclibpy.media import ThermodynamicState + +logger = logging.getLogger(__name__) + + +def calc_area(Q, k, Tm): + if Tm * k == 0: + return 1e10 # A really large number, but not np.inf to still be able to calculate with it + else: + return Q / (k * Tm) + + +def calc_mean_temperature(T_hot_in, T_hot_out, T_cold_in, T_cold_out): + dT_A = T_hot_in - T_cold_out + dT_B = T_hot_out - T_cold_in + if dT_B < 0 or dT_A < 0: + return 0 # Heat can't be transferred + if np.isclose(dT_B, 0) or np.isclose(dT_A, 0) or np.isclose(dT_B, dT_A): + return abs((dT_A + dT_B) / 2) + return (dT_A - dT_B) / np.log(dT_A / dT_B) + + +def separate_phases(m_flow, med_prop, state_max: ThermodynamicState, state_min: ThermodynamicState, p: float): + """ + Separates a flow with possible phase changes into three parts: + subcooling (sc), latent phase change (lat), and superheating (sh) + at the given pressure. + + Args: + state_max (ThermodynamicState): State with higher enthalpy. + state_min (ThermodynamicState): State with lower enthalpy. + p (float): Pressure of phase change. + + Returns: + Tuple[float, float, float, ThermodynamicState, ThermodynamicState]: + Q_sc: Heat for subcooling. + Q_lat: Heat for latent phase change. + Q_sh: Heat for superheating. + state_q0: State at vapor quality 0 and the given pressure. + state_q1: State at vapor quality 1 and the given pressure. + """ + # Get relevant states: + state_q0 = med_prop.calc_state("PQ", p, 0) + state_q1 = med_prop.calc_state("PQ", p, 1) + Q_sc = max(0.0, + min((state_q0.h - state_min.h), + (state_max.h - state_min.h))) * m_flow + Q_lat = max(0.0, + (min(state_max.h, state_q1.h) - + max(state_min.h, state_q0.h))) * m_flow + Q_sh = max(0.0, + min((state_max.h - state_q1.h), + (state_max.h - state_min.h))) * m_flow + return Q_sc, Q_lat, Q_sh, state_q0, state_q1 + + +class MovingBoundaryTmCondenser(HeatExchanger): + """ + Condenser class which implements the actual `calc` method. + + Assumptions: + - No phase changes in secondary medium + - cp of secondary medium is constant over heat-exchanger + + See parent classes for arguments. + """ + + def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): + """ + Calculate the heat exchanger with the Tm-Method based on the given inputs. + + The flowsheet state can be used to save important variables + during calculation for later analysis. + + Both return values are used to check if the heat transfer is valid or not. + + Args: + inputs (Inputs): The inputs for the calculation. + fs_state (FlowsheetState): The flowsheet state to save important variables. + + Returns: + Tuple[float, float]: + error: Error in percentage between the required and calculated heat flow rates. + dT_min: Minimal temperature difference (can be negative). + """ + self.m_flow_secondary = inputs.m_flow_con # [kg/s] + self.calc_secondary_cp(T=inputs.T_con_in) + + # First we separate the flow: + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_inlet, + state_min=self.state_outlet, + p=self.state_inlet.p + ) + Q = Q_sc + Q_lat + Q_sh + + # Note: As Q_con_Tm has to converge to Q_con (m_ref*delta_h), we can safely + # calculate the output temperature. + + T_mean = inputs.T_con_in + self.calc_secondary_Q_flow(Q) / (self.m_flow_secondary_cp * 2) + tra_prop_med = self.calc_transport_properties_secondary_medium(T_mean) + alpha_med_wall = self.calc_alpha_secondary(tra_prop_med) + + # Calculate secondary_medium side temperatures: + # Assumption loss is the same correlation for each regime + T_sc = inputs.T_con_in + self.calc_secondary_Q_flow(Q_sc) / self.m_flow_secondary_cp + T_sh = T_sc + self.calc_secondary_Q_flow(Q_lat) / self.m_flow_secondary_cp + T_out = T_sh + self.calc_secondary_Q_flow(Q_sh) / self.m_flow_secondary_cp + + # 1. Regime: Subcooling + Q_sc_Tm, A_sc, A_sc_available = 0, 0, 0 + if Q_sc > 0 and (state_q0.T != self.state_outlet.T): + # Get transport properties: + tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(state_q0, self.state_outlet) + alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_con) + + # Only use still available area: + k_sc = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sc = calc_mean_temperature( + T_hot_in=state_q0.T, T_hot_out=self.state_outlet.T, + T_cold_in=inputs.T_con_in, T_cold_out=T_sc + ) + A_sc = calc_area(Q_sc, k_sc, T_m_sc) + A_sc_available = min(self.A, A_sc) + Q_sc_Tm = A_sc_available * k_sc * T_m_sc + + # 2. Regime: Latent heat exchange + Q_lat_Tm, A_lat, A_lat_available = 0, 0, 0 + if Q_lat > 0: + # Get transport properties: + alpha_ref_wall = self.calc_alpha_two_phase( + state_q0=state_q0, + state_q1=state_q1, + fs_state=fs_state, + inputs=inputs + ) + k_lat = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_lat = calc_mean_temperature( + T_hot_in=state_q1.T, T_hot_out=state_q0.T, + T_cold_in=T_sc, T_cold_out=T_sh + ) + A_lat = calc_area(Q_lat, k_lat, T_m_lat) + # Only use still available area: + A_lat_available = min(max(self.A - A_sc_available, 0), A_lat) + Q_lat_Tm = A_lat_available * k_lat * T_m_lat + + logger.debug(f"con_lat: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + # 3. Regime: Superheat heat exchange + Q_sh_Tm, A_sh = 0, 0 + if Q_sh and (self.state_inlet.T != state_q1.T): + # Get transport properties: + tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(self.state_inlet, state_q1) + alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_con) + + k_sh = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sh = calc_mean_temperature( + T_hot_in=self.state_inlet.T, T_hot_out=state_q1.T, + T_cold_in=T_sh, T_cold_out=T_out + ) + A_sh = calc_area(Q_sh, k_sh, T_m_sh) + # Only use still available area: + A_sh_available = min(max(self.A - A_sc_available - A_lat_available, 0), A_sh) + Q_sh_Tm = A_sh_available * k_sh * T_m_sh + + logger.debug(f"con_sh: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + A_necessary = A_sh + A_lat + A_sc + Q_Tm = Q_sh_Tm + Q_sc_Tm + Q_lat_Tm + error_A = (1 - A_necessary / self.A) * 100 + error = (Q_Tm / Q - 1) * 100 + # Get possible dT_min: + dT_min_in = self.state_outlet.T - inputs.T_con_in + dT_min_out = self.state_inlet.T - T_out + dT_min_LatSH = state_q1.T - T_sh + + fs_state.set(name="A_con_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in condenser") + fs_state.set(name="A_con_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in condenser") + fs_state.set(name="A_con_sc", value=A_sc, unit="m2", description="Area for subcooling heat exchange in condenser") + + return error, min(dT_min_in, + dT_min_LatSH, + dT_min_out) + + +class MovingBoundaryTmEvaporator(HeatExchanger): + """ + Evaporator class which implements the actual `calc` method. + + Assumptions: + - No phase changes in secondary medium + - cp of secondary medium is constant over heat-exchanger + + See parent classes for arguments. + """ + + def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): + """ + Calculate the heat exchanger with the Tm-Method based on the given inputs. + + The flowsheet state can be used to save important variables + during calculation for later analysis. + + Both return values are used to check if the heat transfer is valid or not. + + Args: + inputs (Inputs): The inputs for the calculation. + fs_state (FlowsheetState): The flowsheet state to save important variables. + + Returns: + Tuple[float, float]: + error: Error in percentage between the required and calculated heat flow rates. + dT_min: Minimal temperature difference (can be negative). + """ + self.m_flow_secondary = inputs.m_flow_eva # [kg/s] + self.calc_secondary_cp(T=inputs.T_eva_in) + + # First we separate the flow: + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_outlet, + state_min=self.state_inlet, + p=self.state_outlet.p + ) + + Q = Q_sc + Q_lat + Q_sh + + # Note: As Q_eva_Tm has to converge to Q_eva (m_ref*delta_h), we can safely + # calculate the output temperature. + T_mean = inputs.T_eva_in - Q / (self.m_flow_secondary_cp * 2) + tra_prop_med = self.calc_transport_properties_secondary_medium(T_mean) + alpha_med_wall = self.calc_alpha_secondary(tra_prop_med) + # alpha_med_wall = 26 + + # Calculate secondary_medium side temperatures: + T_sh = inputs.T_eva_in - Q_sh / self.m_flow_secondary_cp + T_sc = T_sh - Q_lat / self.m_flow_secondary_cp + T_out = T_sc - Q_sc / self.m_flow_secondary_cp + + # 1. Regime: Superheating + Q_sh_Tm, A_sh, A_sh_available = 0, 0, 0 + if Q_sh and (self.state_outlet.T != state_q1.T): + # Get transport properties: + tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(self.state_outlet, state_q1) + alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_eva) + + k_sh = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sh = calc_mean_temperature( + T_hot_in=inputs.T_eva_in, T_hot_out=T_sh, + T_cold_in=state_q1.T, T_cold_out=self.state_outlet.T + ) + # Only use still available area: + A_sh = calc_area(Q_sh, k_sh, T_m_sh) + A_sh_available = min(self.A, A_sh) + Q_sh_Tm = A_sh_available * k_sh * T_m_sh + + logger.debug(f"eva_sh: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + # 2. Regime: Latent heat exchange + Q_lat_Tm, A_lat, A_lat_available = 0, 0, 0 + if Q_lat > 0: + alpha_ref_wall = self.calc_alpha_two_phase( + state_q0=state_q0, + state_q1=state_q1, + fs_state=fs_state, + inputs=inputs + ) + k_lat = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_lat = calc_mean_temperature( + T_hot_in=T_sh, T_hot_out=T_sc, + T_cold_in=state_q0.T, T_cold_out=state_q1.T + ) + A_lat = calc_area(Q_lat, k_lat, T_m_lat) + # Only use still available area: + A_lat_available = min(max(self.A - A_sh_available, 0), A_lat) + Q_lat_Tm = A_lat_available * k_lat * T_m_lat + + logger.debug(f"eva_lat: pri: {round(alpha_ref_wall, 2)} sec: {round(alpha_med_wall, 2)}") + + # 3. Regime: Subcooling + Q_sc_Tm, A_sc = 0, 0 + if Q_sc > 0 and (state_q0.T != self.state_inlet.T): + # Get transport properties: + tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(state_q0, self.state_inlet) + alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_eva) + + # Only use still available area: + k_sc = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) + T_m_sc = calc_mean_temperature( + T_hot_in=T_sc, T_hot_out=T_out, + T_cold_in=self.state_inlet.T, T_cold_out=state_q0.T + ) + A_sc = calc_area(Q_sc, k_sc, T_m_sc) + A_sc_available = min(max(self.A - A_sh_available - A_lat_available, 0), A_sc) + Q_sc_Tm = A_sc_available * k_sc * T_m_sc + + A_necessary = A_sc + A_lat + A_sh + Q_Tm = Q_sh_Tm + Q_sc_Tm + Q_lat_Tm + error_A = (A_necessary / self.A - 1) * 100 + error = (Q_Tm / Q - 1) * 100 + # Get dT_min + dT_min_in = inputs.T_eva_in - self.state_outlet.T + dT_min_out = T_out - self.state_inlet.T + + fs_state.set(name="A_eva_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in evaporator") + fs_state.set(name="A_eva_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in evaporator") + + return error, min(dT_min_out, dT_min_in) diff --git a/vclibpy/components/heat_exchangers/ntu.py b/vclibpy/components/heat_exchangers/ntu.py index db7c92b..e1fe9de 100644 --- a/vclibpy/components/heat_exchangers/ntu.py +++ b/vclibpy/components/heat_exchangers/ntu.py @@ -28,7 +28,7 @@ class BasicNTU(HeatExchanger, abc.ABC): Typical values are around 20-30. """ - def __init__(self, flow_type: str, ratio_outer_to_inner_area: float, **kwargs): + def __init__(self, flow_type: str, **kwargs): """ Initializes BasicNTU. @@ -45,7 +45,6 @@ def __init__(self, flow_type: str, ratio_outer_to_inner_area: float, **kwargs): **kwargs: Additional keyword arguments passed to the parent class. """ super().__init__(**kwargs) - self.ratio_outer_to_inner_area = ratio_outer_to_inner_area # Set primary cp: self._primary_cp = None @@ -102,26 +101,6 @@ def calc_R(self) -> float: return self.m_flow_secondary_cp / m_flow_pri_cp return m_flow_pri_cp / self.m_flow_secondary_cp - def calc_k(self, alpha_pri: float, alpha_sec: float) -> float: - """ - Calculate the overall heat transfer coefficient (k) of the heat exchanger. - - Args: - alpha_pri (float): Heat transfer coefficient for the primary medium. - alpha_sec (float): Heat transfer coefficient for the secondary medium. - - Returns: - float: Overall heat transfer coefficient (k). - """ - k_wall = self.calc_wall_heat_transfer() - k = (1 / ( - (1 / alpha_pri) * self.ratio_outer_to_inner_area + - (1 / k_wall) * self.ratio_outer_to_inner_area + - (1 / alpha_sec) - ) - ) - return k - @staticmethod def calc_NTU(A: float, k: float, m_flow_cp: float) -> float: """ From 92c85308f0e6c5b6566d982be7d909a35f2dbebf Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 13:48:53 +0100 Subject: [PATCH 05/36] Separate functions, improve iteration analysis plotting --- vclibpy/flowsheets/base.py | 99 +++++++++++++++++++++++++++++--------- 1 file changed, 77 insertions(+), 22 deletions(-) diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index b01fca2..2541b2f 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -71,6 +71,15 @@ def terminate(self): def get_all_components(self) -> List[BaseComponent]: return [self.condenser, self.evaporator] + def get_start_condensing_pressure(self, inputs: Inputs): + T_3_start = inputs.T_con_in + inputs.dT_con_subcooling + p_2_start = self.med_prop.calc_state("TQ", T_3_start, 0).p + return p_2_start + + def get_start_evaporating_pressure(self, inputs: Inputs, dT_pinch: float = 0): + T_1_start = inputs.T_eva_in - inputs.dT_eva_superheating - dT_pinch + return self.med_prop.calc_state("TQ", T_1_start, 1).p + def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): """ Calculate the steady-state performance of a vapor compression cycle @@ -126,15 +135,21 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): # Settings min_iteration_step = kwargs.pop("min_iteration_step", 1) save_path_plots = kwargs.get("save_path_plots", None) - input_name = ";".join([k + "=" + str(np.round(v.value, 3)).replace(".", "_") - for k, v in inputs.get_variables().items()]) show_iteration = kwargs.get("show_iteration", False) use_quick_solver = kwargs.pop("use_quick_solver", True) err_ntu = kwargs.pop("max_err_ntu", 0.5) err_dT_min = kwargs.pop("max_err_dT_min", 0.1) max_num_iterations = kwargs.pop("max_num_iterations", 1e5) + dT_pinch_con_guess = kwargs.pop("dT_pinch_con_guess", 0) + dT_pinch_eva_guess = kwargs.pop("dT_pinch_eva_guess", 0) p_1_history = [] p_2_history = [] + error_con_history = [] + error_eva_history = [] + dT_eva_history = [] + dT_con_history = [] + error_con, dT_min_eva, dT_min_con, error_eva = np.NAN, np.NAN, np.NAN, np.NAN + plot_last = -100 if use_quick_solver: step_p1 = kwargs.get("step_max", 10000) @@ -149,19 +164,15 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): self.setup_new_fluid(fluid) # First: Iterate with given conditions to get the 4 states and the mass flow rate: - T_1_start = inputs.T_eva_in - inputs.dT_eva_superheating - T_3_start = inputs.T_con_in + inputs.dT_con_subcooling - p_1_start = self.med_prop.calc_state("TQ", T_1_start, 1).p - p_2_start = self.med_prop.calc_state("TQ", T_3_start, 0).p - p_1_next = p_1_start - p_2_next = p_2_start + p_1_next = self.get_start_evaporating_pressure(inputs=inputs, dT_pinch=dT_pinch_eva_guess) + p_2_next = self.get_start_condensing_pressure(inputs=inputs) fs_state = FlowsheetState() # Always log what is happening in the whole flowsheet fs_state.set(name="Q_con", value=1, unit="W", description="Condenser heat flow rate") fs_state.set(name="COP", value=0, unit="-", description="Coefficient of performance") if show_iteration: - fig_iterations, ax_iterations = plt.subplots(2) + fig_iterations, ax_iterations = plt.subplots(3, 2, sharex=True) num_iterations = 0 @@ -178,13 +189,29 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): p_1 = p_1_next p_2 = p_2_next - p_1_history.append(p_1) - p_2_history.append(p_2) + p_1_history.append(p_1 / 1e5) + p_2_history.append(p_2 / 1e5) + error_con_history.append(error_con) + error_eva_history.append(error_eva) + dT_con_history.append(dT_min_con) + dT_eva_history.append(dT_min_eva) + if show_iteration: - ax_iterations[0].cla() - ax_iterations[1].cla() - ax_iterations[0].scatter(list(range(len(p_1_history))), p_1_history) - ax_iterations[1].scatter(list(range(len(p_2_history))), p_2_history) + for ax in ax_iterations.flatten(): + ax.clear() + iterations = list(range(len(p_1_history)))[plot_last:] + ax_iterations[0, 0].set_ylabel("error_eva in %") + ax_iterations[0, 1].set_ylabel("error_con in %") + ax_iterations[1, 0].set_ylabel("$\Delta T_\mathrm{Min}$ in K") + ax_iterations[1, 1].set_ylabel("$\Delta T_\mathrm{Min}$ in K") + ax_iterations[2, 0].set_ylabel("$p_1$ in bar") + ax_iterations[2, 1].set_ylabel("$p_2$ in bar") + ax_iterations[0, 0].scatter(iterations, error_eva_history[plot_last:]) + ax_iterations[0, 1].scatter(iterations, error_con_history[plot_last:]) + ax_iterations[1, 0].scatter(iterations, dT_eva_history[plot_last:]) + ax_iterations[1, 1].scatter(iterations, dT_con_history[plot_last:]) + ax_iterations[2, 0].scatter(iterations, p_1_history[plot_last:]) + ax_iterations[2, 1].scatter(iterations, p_2_history[plot_last:]) plt.draw() plt.pause(1e-5) @@ -206,18 +233,24 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): step_p1 /= 10 continue - # Calculate the states based on the given flowsheet try: - self.calc_states(p_1, p_2, inputs=inputs, fs_state=fs_state) + error_eva, dT_min_eva, error_con, dT_min_con = self.calculate_cycle_for_pressures( + p_1=p_1, p_2=p_2, inputs=inputs, fs_state=fs_state + ) except ValueError as err: logger.error("An error occurred while calculating states. " "Can't guess next pressures, thus, exiting: %s", err) return - if save_path_plots is not None and num_iterations == 1 and show_iteration: - self.plot_cycle(save_path=save_path_plots.joinpath(f"{input_name}_initialization.png"), inputs=inputs) - # Check heat exchangers: - error_eva, dT_min_eva = self.evaporator.calc(inputs=inputs, fs_state=fs_state) + if num_iterations == 1: + + if save_path_plots is not None and show_iteration: + input_name = inputs.get_name() + self.plot_cycle( + save_path=save_path_plots.joinpath(f"{input_name}_initialization.png"), + inputs=inputs + ) + if not isinstance(error_eva, float): print(error_eva) if error_eva < 0: @@ -233,7 +266,6 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): p_1_next = p_1 + step_p1 continue - error_con, dT_min_con = self.condenser.calc(inputs=inputs, fs_state=fs_state) if error_con < 0: p_2_next = p_2 + step_p2 continue @@ -273,6 +305,28 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): if show_iteration: plt.close(fig_iterations) + return self.calculate_outputs_for_valid_pressures( + p_1=p_1, p_2=p_2, inputs=inputs, fs_state=fs_state, + save_path_plots=save_path_plots + ) + + def calculate_cycle_for_pressures(self, p_1: float, p_2: float, inputs: Inputs, fs_state: FlowsheetState): + # Calculate the states based on the given flowsheet + self.calc_states(p_1, p_2, inputs=inputs, fs_state=fs_state) + # Check heat exchangers: + error_eva, dT_min_eva = self.evaporator.calc(inputs=inputs, fs_state=fs_state) + error_con, dT_min_con = self.condenser.calc(inputs=inputs, fs_state=fs_state) + return error_eva, dT_min_eva, error_con, dT_min_con + + def calculate_outputs_for_valid_pressures( + self, + p_1, + p_2, + fs_state: FlowsheetState, + inputs: Inputs, + save_path_plots + ): + self.calc_states(p_1, p_2, inputs=inputs, fs_state=fs_state) # Calculate the heat flow rates for the selected states. Q_con = self.condenser.calc_Q_flow() Q_con_outer = self.condenser.calc_secondary_Q_flow(Q_con) @@ -321,6 +375,7 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): ) if save_path_plots is not None: + input_name = inputs.get_name() self.plot_cycle(save_path=save_path_plots.joinpath(f"{input_name}_final_result.png"), inputs=inputs) return fs_state From 3ba9f65a5c5fbd4d9878eafa44a4248775d7f8bc Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 13:51:15 +0100 Subject: [PATCH 06/36] Add fsolve based solvers and first iteration guess for condenser --- vclibpy/flowsheets/base.py | 86 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index 2541b2f..f607937 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -1,6 +1,7 @@ import logging from typing import List import numpy as np +from scipy.optimize import fsolve from abc import abstractmethod import matplotlib.pyplot as plt @@ -76,6 +77,26 @@ def get_start_condensing_pressure(self, inputs: Inputs): p_2_start = self.med_prop.calc_state("TQ", T_3_start, 0).p return p_2_start + def improve_first_condensing_guess(self, inputs: Inputs, m_flow_guess, p_2_guess, dT_pinch_assumption=0): + self.condenser.m_flow_secondary = inputs.m_flow_con # [kg/s] + self.condenser.calc_secondary_cp(T=inputs.T_con_in) + + def nonlinear_func(p_2, *args): + _flowsheet, _inputs, _m_flow_ref, _dT_pinch = args + state_q0 = _flowsheet.med_prop.calc_state("PQ", p_2, 0) + state_q1 = _flowsheet.med_prop.calc_state("PQ", p_2, 1) + state_3 = _flowsheet.med_prop.calc_state("PT", p_2, state_q0.T - _inputs.dT_con_subcooling) + Q_water_till_q1 = (state_q1.h - state_3.h) * _m_flow_ref + T_water_q1 = _inputs.T_con_in + Q_water_till_q1 / _flowsheet.condenser.m_flow_secondary_cp + return T_water_q1 + _dT_pinch - state_q1.T + + p_2_guess_optimized = fsolve( + func=nonlinear_func, + x0=p_2_guess, + args=(self, inputs, m_flow_guess, dT_pinch_assumption) + )[0] + return p_2_guess_optimized + def get_start_evaporating_pressure(self, inputs: Inputs, dT_pinch: float = 0): T_1_start = inputs.T_eva_in - inputs.dT_eva_superheating - dT_pinch return self.med_prop.calc_state("TQ", T_1_start, 1).p @@ -142,6 +163,7 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): max_num_iterations = kwargs.pop("max_num_iterations", 1e5) dT_pinch_con_guess = kwargs.pop("dT_pinch_con_guess", 0) dT_pinch_eva_guess = kwargs.pop("dT_pinch_eva_guess", 0) + improve_first_condensing_guess = kwargs.pop("improve_first_condensing_guess", False) p_1_history = [] p_2_history = [] error_con_history = [] @@ -243,6 +265,13 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): return if num_iterations == 1: + if improve_first_condensing_guess: + p_2_next = self.improve_first_condensing_guess( + inputs=inputs, + m_flow_guess=self.condenser.m_flow, + p_2_guess=p_2, + dT_pinch_assumption=dT_pinch_con_guess + ) if save_path_plots is not None and show_iteration: input_name = inputs.get_name() @@ -310,6 +339,63 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): save_path_plots=save_path_plots ) + def calc_steady_state_fsolve(self, inputs: Inputs, fluid: str = None, **kwargs): + # Settings + max_err = kwargs.pop("max_err_ntu", 0.5) + + save_path_plots = kwargs.get("save_path_plots", None) + dT_pinch_eva_guess = kwargs.pop("dT_pinch_eva_guess", 0) + + # Setup fluid: + if fluid is None: + fluid = self.fluid + self.setup_new_fluid(fluid) + + # First: Iterate with given conditions to get the 4 states and the mass flow rate: + p_1_next = self.get_start_evaporating_pressure(inputs=inputs, dT_pinch=dT_pinch_eva_guess) + p_2_next = self.get_start_condensing_pressure(inputs=inputs) + + def nonlinear_func(x, *args): + _flowsheet, _inputs, _fs_state, _max_err = args + _p_1, _p_2 = x + if not (_p_1 < _p_2 < self._p_max): + return 1000, 1000 + if not (self._p_min < _p_1 < _p_2): + return 1000, 1000 + _error_eva, dT_min_eva, _error_con, dT_min_con = _flowsheet.calculate_cycle_for_pressures( + p_1=x[0], p_2=x[1], inputs=_inputs, fs_state=_fs_state + ) + + if 0 <= _error_eva < _max_err: + _error_eva = 0 + if 0 <= _error_con < _max_err: + _error_con = 0 + if 0 > _error_eva: + _error_eva *= 5 + if 0 > _error_con: + _error_con *= 5 + return _error_eva, _error_con + + fs_state = FlowsheetState() # Always log what is happening in the whole flowsheet + try: + args = (self, inputs, fs_state, max_err) + p_optimized = fsolve( + func=nonlinear_func, + x0=np.array([p_1_next, p_2_next]), + args=args + ) + except Exception as err: + logger.error("An error occurred while calculating states using fsolve: %s", err) + return + error_con, error_eva = nonlinear_func(p_optimized, *args) + print(f"{error_con=}, {error_eva=}") + + p_1, p_2 = p_optimized + return self.calculate_outputs_for_valid_pressures( + p_1=p_1, p_2=p_2, inputs=inputs, fs_state=fs_state, + save_path_plots=save_path_plots + ) + def calculate_cycle_for_pressures(self, p_1: float, p_2: float, inputs: Inputs, fs_state: FlowsheetState): # Calculate the states based on the given flowsheet self.calc_states(p_1, p_2, inputs=inputs, fs_state=fs_state) From 65118d6b17673e3eb804a62fed318e88b7b2f384 Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 13:59:06 +0100 Subject: [PATCH 07/36] add errors to fs_state --- vclibpy/flowsheets/base.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index f607937..79621ea 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -418,8 +418,8 @@ def calculate_outputs_for_valid_pressures( Q_con_outer = self.condenser.calc_secondary_Q_flow(Q_con) Q_eva = self.evaporator.calc_Q_flow() Q_eva_outer = self.evaporator.calc_secondary_Q_flow(Q_eva) - self.evaporator.calc(inputs=inputs, fs_state=fs_state) - self.condenser.calc(inputs=inputs, fs_state=fs_state) + error_con, dT_min_con = self.condenser.calc(inputs=inputs, fs_state=fs_state) + error_eva, dT_min_eva = self.evaporator.calc(inputs=inputs, fs_state=fs_state) P_el = self.calc_electrical_power(fs_state=fs_state, inputs=inputs) T_con_out = inputs.T_con_in + Q_con_outer / self.condenser.m_flow_secondary_cp @@ -459,7 +459,22 @@ def calculate_outputs_for_valid_pressures( name="COP_outer", value=COP_outer, unit="-", description="Outer COP, including heat losses" ) - + fs_state.set( + name="error_con", value=error_con, + unit="%", description="Error in condenser heat exchanger model" + ) + fs_state.set( + name="error_eva", value=error_eva, + unit="%", description="Error in evaporator heat exchanger model" + ) + fs_state.set( + name="dT_min_eva", value=dT_min_eva, + unit="K", description="Evaporator pinch temperature" + ) + fs_state.set( + name="dT_min_con", value=dT_min_con, + unit="K", description="Condenser pinch temperature" + ) if save_path_plots is not None: input_name = inputs.get_name() self.plot_cycle(save_path=save_path_plots.joinpath(f"{input_name}_final_result.png"), inputs=inputs) From 621a4b450710ecffacafcf35a0a5b8bd0743e8ca Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Thu, 21 Dec 2023 21:07:37 +0100 Subject: [PATCH 08/36] add flow_type to heat exchanger --- .../components/heat_exchangers/heat_exchanger.py | 1 + .../heat_exchangers/moving_boundary_Tm.py | 16 +++++++++++++--- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/vclibpy/components/heat_exchangers/heat_exchanger.py b/vclibpy/components/heat_exchangers/heat_exchanger.py index 579ed49..ebb8f5e 100644 --- a/vclibpy/components/heat_exchangers/heat_exchanger.py +++ b/vclibpy/components/heat_exchangers/heat_exchanger.py @@ -39,6 +39,7 @@ def __init__( two_phase_heat_transfer: TwoPhaseHeatTransfer, secondary_medium: str, ratio_outer_to_inner_area: float = 1, + flow_type: str = "counter" ): super().__init__() self.A = A diff --git a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py index 1ebec10..76265cf 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py @@ -15,9 +15,19 @@ def calc_area(Q, k, Tm): return Q / (k * Tm) -def calc_mean_temperature(T_hot_in, T_hot_out, T_cold_in, T_cold_out): - dT_A = T_hot_in - T_cold_out - dT_B = T_hot_out - T_cold_in +def calc_mean_temperature( + T_hot_in: float, T_hot_out: float, + T_cold_in: float, T_cold_out: float, + flow_type: str = "counter" +): + if flow_type == "counter": + dT_A = T_hot_in - T_cold_out + dT_B = T_hot_out - T_cold_in + elif flow_type == "parallel": + dT_A = T_hot_in - T_cold_in + dT_B = T_hot_out - T_cold_out + else: + raise TypeError("Given flow_type is not supported yet") if dT_B < 0 or dT_A < 0: return 0 # Heat can't be transferred if np.isclose(dT_B, 0) or np.isclose(dT_A, 0) or np.isclose(dT_B, dT_A): From 42aba5d3eda3c60064bfe71f8d2858d11c239b84 Mon Sep 17 00:00:00 2001 From: FWuellhorst Date: Wed, 28 Feb 2024 17:53:42 +0100 Subject: [PATCH 09/36] fix str name --- vclibpy/components/heat_exchangers/heat_exchanger.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vclibpy/components/heat_exchangers/heat_exchanger.py b/vclibpy/components/heat_exchangers/heat_exchanger.py index ebb8f5e..23e8fbe 100644 --- a/vclibpy/components/heat_exchangers/heat_exchanger.py +++ b/vclibpy/components/heat_exchangers/heat_exchanger.py @@ -64,7 +64,7 @@ def start_secondary_med_prop(self): # Set up the secondary_medium wrapper: med_prop_class, med_prop_kwargs = media.get_global_med_prop_and_kwargs() if self.secondary_medium == "air" and med_prop_class == media.RefProp: - fluid_name = "AIR.PPF" + fluid_name = "air.ppf" else: fluid_name = self.secondary_medium if self.med_prop_sec is not None: From 8fd7bc8f70dfab7bd9c8bcad792db5f5cbd1471c Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Wed, 15 May 2024 12:05:53 +0200 Subject: [PATCH 10/36] add option for T_con_out usage in all simulations --- examples/e6_simple_heat_pump.py | 7 +- examples/e7_vapor_injection.py | 7 +- tests/test_regression.py | 4 +- .../heat_exchangers/moving_boundary_Tm.py | 113 ++------- .../heat_exchangers/moving_boundary_ntu.py | 230 ++++++++---------- vclibpy/components/heat_exchangers/utils.py | 93 +++++++ vclibpy/datamodels.py | 39 ++- vclibpy/flowsheets/base.py | 45 +++- vclibpy/utils/automation.py | 84 +++++-- 9 files changed, 358 insertions(+), 264 deletions(-) create mode 100644 vclibpy/components/heat_exchangers/utils.py diff --git a/examples/e6_simple_heat_pump.py b/examples/e6_simple_heat_pump.py index e258898..7714f7e 100644 --- a/examples/e6_simple_heat_pump.py +++ b/examples/e6_simple_heat_pump.py @@ -61,9 +61,11 @@ def main(): ) # As in the other example, we can specify save-paths, # solver settings and inputs to vary: + # Note that T_con can either be inlet or outlet, depending on the setting + # of `use_condenser_inlet` later on. In this case, we simulate the inlet, T_con_in save_path = r"D:\00_temp\simple_heat_pump" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] - T_con_in_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] + T_con_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] n_ar = [0.3, 0.7, 1] # Now, we can generate the full-factorial performance map @@ -77,9 +79,10 @@ def main(): save_path_sdf, save_path_csv = utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=True, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, diff --git a/examples/e7_vapor_injection.py b/examples/e7_vapor_injection.py index 315b801..ca2a811 100644 --- a/examples/e7_vapor_injection.py +++ b/examples/e7_vapor_injection.py @@ -69,7 +69,7 @@ def main(): # solver settings and inputs to vary: save_path = r"D:\00_temp\vapor_injection" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] - T_con_in_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] + T_con_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] n_ar = [0.3, 0.7, 1] # Now, we can generate the full-factorial performance map @@ -83,9 +83,10 @@ def main(): utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=True, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -127,7 +128,7 @@ def main(): utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=r"D:\00_temp\vapor_injection_economizer", - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, use_multiprocessing=False, diff --git a/tests/test_regression.py b/tests/test_regression.py index 90196c5..3501e0e 100644 --- a/tests/test_regression.py +++ b/tests/test_regression.py @@ -131,7 +131,7 @@ def _regression_of_examples(self, flowsheet, fluid): # Just for quick study: Specify concrete points: T_eva_in_ar = [-10 + 273.15, 273.15] - T_con_in_ar = [30 + 273.15, 70 + 273.15] + T_con_ar = [30 + 273.15, 70 + 273.15] n_ar = [0.3, 1] os.makedirs(self.working_dir, exist_ok=True) @@ -146,7 +146,7 @@ def _regression_of_examples(self, flowsheet, fluid): _, path_csv = utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=self.working_dir, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, use_multiprocessing=False, diff --git a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py index 76265cf..649d432 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py @@ -1,74 +1,12 @@ import logging -import numpy as np +from vclibpy.components.heat_exchangers import utils from vclibpy.datamodels import FlowsheetState, Inputs from vclibpy.components.heat_exchangers import HeatExchanger -from vclibpy.media import ThermodynamicState logger = logging.getLogger(__name__) -def calc_area(Q, k, Tm): - if Tm * k == 0: - return 1e10 # A really large number, but not np.inf to still be able to calculate with it - else: - return Q / (k * Tm) - - -def calc_mean_temperature( - T_hot_in: float, T_hot_out: float, - T_cold_in: float, T_cold_out: float, - flow_type: str = "counter" -): - if flow_type == "counter": - dT_A = T_hot_in - T_cold_out - dT_B = T_hot_out - T_cold_in - elif flow_type == "parallel": - dT_A = T_hot_in - T_cold_in - dT_B = T_hot_out - T_cold_out - else: - raise TypeError("Given flow_type is not supported yet") - if dT_B < 0 or dT_A < 0: - return 0 # Heat can't be transferred - if np.isclose(dT_B, 0) or np.isclose(dT_A, 0) or np.isclose(dT_B, dT_A): - return abs((dT_A + dT_B) / 2) - return (dT_A - dT_B) / np.log(dT_A / dT_B) - - -def separate_phases(m_flow, med_prop, state_max: ThermodynamicState, state_min: ThermodynamicState, p: float): - """ - Separates a flow with possible phase changes into three parts: - subcooling (sc), latent phase change (lat), and superheating (sh) - at the given pressure. - - Args: - state_max (ThermodynamicState): State with higher enthalpy. - state_min (ThermodynamicState): State with lower enthalpy. - p (float): Pressure of phase change. - - Returns: - Tuple[float, float, float, ThermodynamicState, ThermodynamicState]: - Q_sc: Heat for subcooling. - Q_lat: Heat for latent phase change. - Q_sh: Heat for superheating. - state_q0: State at vapor quality 0 and the given pressure. - state_q1: State at vapor quality 1 and the given pressure. - """ - # Get relevant states: - state_q0 = med_prop.calc_state("PQ", p, 0) - state_q1 = med_prop.calc_state("PQ", p, 1) - Q_sc = max(0.0, - min((state_q0.h - state_min.h), - (state_max.h - state_min.h))) * m_flow - Q_lat = max(0.0, - (min(state_max.h, state_q1.h) - - max(state_min.h, state_q0.h))) * m_flow - Q_sh = max(0.0, - min((state_max.h - state_q1.h), - (state_max.h - state_min.h))) * m_flow - return Q_sc, Q_lat, Q_sh, state_q0, state_q1 - - class MovingBoundaryTmCondenser(HeatExchanger): """ Condenser class which implements the actual `calc` method. @@ -99,10 +37,10 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): dT_min: Minimal temperature difference (can be negative). """ self.m_flow_secondary = inputs.m_flow_con # [kg/s] - self.calc_secondary_cp(T=inputs.T_con_in) + self.calc_secondary_cp(T=inputs.T_con) # First we separate the flow: - Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = utils.separate_phases( m_flow=self.m_flow, med_prop=self.med_prop, state_max=self.state_inlet, @@ -111,19 +49,14 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): ) Q = Q_sc + Q_lat + Q_sh - # Note: As Q_con_Tm has to converge to Q_con (m_ref*delta_h), we can safely - # calculate the output temperature. + T_in, T_sc, T_sh, T_out = utils.get_condenser_phase_temperatures_and_alpha( + heat_exchanger=self, inputs=inputs, + Q_sc=Q_sc, Q_lat=Q_lat, Q_sh=Q_sh + ) - T_mean = inputs.T_con_in + self.calc_secondary_Q_flow(Q) / (self.m_flow_secondary_cp * 2) - tra_prop_med = self.calc_transport_properties_secondary_medium(T_mean) + tra_prop_med = self.calc_transport_properties_secondary_medium((T_in + T_out) / 2) alpha_med_wall = self.calc_alpha_secondary(tra_prop_med) - # Calculate secondary_medium side temperatures: - # Assumption loss is the same correlation for each regime - T_sc = inputs.T_con_in + self.calc_secondary_Q_flow(Q_sc) / self.m_flow_secondary_cp - T_sh = T_sc + self.calc_secondary_Q_flow(Q_lat) / self.m_flow_secondary_cp - T_out = T_sh + self.calc_secondary_Q_flow(Q_sh) / self.m_flow_secondary_cp - # 1. Regime: Subcooling Q_sc_Tm, A_sc, A_sc_available = 0, 0, 0 if Q_sc > 0 and (state_q0.T != self.state_outlet.T): @@ -133,11 +66,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # Only use still available area: k_sc = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) - T_m_sc = calc_mean_temperature( + T_m_sc = utils.calc_mean_temperature( T_hot_in=state_q0.T, T_hot_out=self.state_outlet.T, - T_cold_in=inputs.T_con_in, T_cold_out=T_sc + T_cold_in=T_in, T_cold_out=T_sc ) - A_sc = calc_area(Q_sc, k_sc, T_m_sc) + A_sc = utils.calc_area(Q_sc, k_sc, T_m_sc) A_sc_available = min(self.A, A_sc) Q_sc_Tm = A_sc_available * k_sc * T_m_sc @@ -152,11 +85,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): inputs=inputs ) k_lat = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) - T_m_lat = calc_mean_temperature( + T_m_lat = utils.calc_mean_temperature( T_hot_in=state_q1.T, T_hot_out=state_q0.T, T_cold_in=T_sc, T_cold_out=T_sh ) - A_lat = calc_area(Q_lat, k_lat, T_m_lat) + A_lat = utils.calc_area(Q_lat, k_lat, T_m_lat) # Only use still available area: A_lat_available = min(max(self.A - A_sc_available, 0), A_lat) Q_lat_Tm = A_lat_available * k_lat * T_m_lat @@ -171,11 +104,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_con) k_sh = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) - T_m_sh = calc_mean_temperature( + T_m_sh = utils.calc_mean_temperature( T_hot_in=self.state_inlet.T, T_hot_out=state_q1.T, T_cold_in=T_sh, T_cold_out=T_out ) - A_sh = calc_area(Q_sh, k_sh, T_m_sh) + A_sh = utils.calc_area(Q_sh, k_sh, T_m_sh) # Only use still available area: A_sh_available = min(max(self.A - A_sc_available - A_lat_available, 0), A_sh) Q_sh_Tm = A_sh_available * k_sh * T_m_sh @@ -187,7 +120,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): error_A = (1 - A_necessary / self.A) * 100 error = (Q_Tm / Q - 1) * 100 # Get possible dT_min: - dT_min_in = self.state_outlet.T - inputs.T_con_in + dT_min_in = self.state_outlet.T - T_in dT_min_out = self.state_inlet.T - T_out dT_min_LatSH = state_q1.T - T_sh @@ -233,7 +166,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): self.calc_secondary_cp(T=inputs.T_eva_in) # First we separate the flow: - Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = utils.separate_phases( m_flow=self.m_flow, med_prop=self.med_prop, state_max=self.state_outlet, @@ -263,12 +196,12 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_eva) k_sh = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) - T_m_sh = calc_mean_temperature( + T_m_sh = utils.calc_mean_temperature( T_hot_in=inputs.T_eva_in, T_hot_out=T_sh, T_cold_in=state_q1.T, T_cold_out=self.state_outlet.T ) # Only use still available area: - A_sh = calc_area(Q_sh, k_sh, T_m_sh) + A_sh = utils.calc_area(Q_sh, k_sh, T_m_sh) A_sh_available = min(self.A, A_sh) Q_sh_Tm = A_sh_available * k_sh * T_m_sh @@ -284,11 +217,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): inputs=inputs ) k_lat = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) - T_m_lat = calc_mean_temperature( + T_m_lat = utils.calc_mean_temperature( T_hot_in=T_sh, T_hot_out=T_sc, T_cold_in=state_q0.T, T_cold_out=state_q1.T ) - A_lat = calc_area(Q_lat, k_lat, T_m_lat) + A_lat = utils.calc_area(Q_lat, k_lat, T_m_lat) # Only use still available area: A_lat_available = min(max(self.A - A_sh_available, 0), A_lat) Q_lat_Tm = A_lat_available * k_lat * T_m_lat @@ -304,11 +237,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # Only use still available area: k_sc = self.calc_k(alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall) - T_m_sc = calc_mean_temperature( + T_m_sc = utils.calc_mean_temperature( T_hot_in=T_sc, T_hot_out=T_out, T_cold_in=self.state_inlet.T, T_cold_out=state_q0.T ) - A_sc = calc_area(Q_sc, k_sc, T_m_sc) + A_sc = utils.calc_area(Q_sc, k_sc, T_m_sc) A_sc_available = min(max(self.A - A_sh_available - A_lat_available, 0), A_sc) Q_sc_Tm = A_sc_available * k_sc * T_m_sc diff --git a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py index d0ee85b..4f68d8b 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py @@ -1,109 +1,69 @@ -import abc import logging import numpy as np from vclibpy.datamodels import FlowsheetState, Inputs from vclibpy.components.heat_exchangers.ntu import BasicNTU -from vclibpy.media import ThermodynamicState +from vclibpy.components.heat_exchangers.utils import separate_phases, get_condenser_phase_temperatures_and_alpha logger = logging.getLogger(__name__) -class MovingBoundaryNTU(BasicNTU, abc.ABC): +def iterate_area(heat_exchanger: BasicNTU, dT_max, alpha_pri, alpha_sec, Q) -> float: """ - Moving boundary NTU based heat exchanger. + Iteratively calculates the required area for the heat exchange. - See parent classe for arguments. - """ - - def separate_phases(self, state_max: ThermodynamicState, state_min: ThermodynamicState, p: float): - """ - Separates a flow with possible phase changes into three parts: - subcooling (sc), latent phase change (lat), and superheating (sh) - at the given pressure. + Args: + heat_exchanger (BasicNTU): An instance of the BasicNTU or children classes + dT_max (float): Maximum temperature differential. + alpha_pri (float): Heat transfer coefficient for the primary medium. + alpha_sec (float): Heat transfer coefficient for the secondary medium. + Q (float): Heat flow rate. - Args: - state_max (ThermodynamicState): State with higher enthalpy. - state_min (ThermodynamicState): State with lower enthalpy. - p (float): Pressure of phase change. - - Returns: - Tuple[float, float, float, ThermodynamicState, ThermodynamicState]: - Q_sc: Heat for subcooling. - Q_lat: Heat for latent phase change. - Q_sh: Heat for superheating. - state_q0: State at vapor quality 0 and the given pressure. - state_q1: State at vapor quality 1 and the given pressure. - """ - # Get relevant states: - state_q0 = self.med_prop.calc_state("PQ", p, 0) - state_q1 = self.med_prop.calc_state("PQ", p, 1) - Q_sc = max(0.0, - min((state_q0.h - state_min.h), - (state_max.h - state_min.h))) * self.m_flow - Q_lat = max(0.0, - (min(state_max.h, state_q1.h) - - max(state_min.h, state_q0.h))) * self.m_flow - Q_sh = max(0.0, - min((state_max.h - state_q1.h), - (state_max.h - state_min.h))) * self.m_flow - return Q_sc, Q_lat, Q_sh, state_q0, state_q1 - - def iterate_area(self, dT_max, alpha_pri, alpha_sec, Q) -> float: - """ - Iteratively calculates the required area for the heat exchange. - - Args: - dT_max (float): Maximum temperature differential. - alpha_pri (float): Heat transfer coefficient for the primary medium. - alpha_sec (float): Heat transfer coefficient for the secondary medium. - Q (float): Heat flow rate. - - Returns: - float: Required area for heat exchange. - """ - _accuracy = 1e-6 # square mm - _step = 1.0 - R = self.calc_R() - k = self.calc_k(alpha_pri, alpha_sec) - m_flow_cp_min = self.calc_m_flow_cp_min() - # First check if point is feasible at all - if dT_max <= 0: - return self.A - eps_necessary = Q / (m_flow_cp_min * dT_max) - - # Special cases: - # --------------- - # eps is equal or higher than 1, an infinite amount of area would be necessary. - if eps_necessary >= 1: - return self.A - # eps is lower or equal to zero: No Area required (Q<=0) - if eps_necessary <= 0: - return 0 - - area = 0.0 - while True: - if self.flow_type == "cross" and area == 0.0: - eps = 0.0 - else: - NTU = self.calc_NTU(area, k, m_flow_cp_min) - eps = self.calc_eps(R, NTU) - if eps >= eps_necessary: - if _step <= _accuracy: - break - else: - # Go back - area -= _step - _step /= 10 - continue - if _step < _accuracy and area > self.A: + Returns: + float: Required area for heat exchange. + """ + _accuracy = 1e-6 # square mm + _step = 1.0 + R = heat_exchanger.calc_R() + k = heat_exchanger.calc_k(alpha_pri, alpha_sec) + m_flow_cp_min = heat_exchanger.calc_m_flow_cp_min() + # First check if point is feasible at all + if dT_max <= 0: + return heat_exchanger.A + eps_necessary = Q / (m_flow_cp_min * dT_max) + + # Special cases: + # --------------- + # eps is equal or higher than 1, an infinite amount of area would be necessary. + if eps_necessary >= 1: + return heat_exchanger.A + # eps is lower or equal to zero: No Area required (Q<=0) + if eps_necessary <= 0: + return 0 + + area = 0.0 + while True: + if heat_exchanger.flow_type == "cross" and area == 0.0: + eps = 0.0 + else: + NTU = heat_exchanger.calc_NTU(area, k, m_flow_cp_min) + eps = heat_exchanger.calc_eps(R, NTU) + if eps >= eps_necessary: + if _step <= _accuracy: break - area += _step + else: + # Go back + area -= _step + _step /= 10 + continue + if _step < _accuracy and area > heat_exchanger.A: + break + area += _step - return min(area, self.A) + return min(area, heat_exchanger.A) -class MovingBoundaryNTUCondenser(MovingBoundaryNTU): +class MovingBoundaryNTUCondenser(BasicNTU): """ Condenser class which implements the actual `calc` method. @@ -133,29 +93,26 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): dT_min: Minimal temperature difference (can be negative). """ self.m_flow_secondary = inputs.m_flow_con # [kg/s] - self.calc_secondary_cp(T=inputs.T_con_in) + self.calc_secondary_cp(T=inputs.T_con) # First we separate the flow: - Q_sc, Q_lat, Q_sh, state_q0, state_q1 = self.separate_phases( - self.state_inlet, - self.state_outlet, - self.state_inlet.p + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_inlet, + state_min=self.state_outlet, + p=self.state_inlet.p ) Q = Q_sc + Q_lat + Q_sh - # Note: As Q_con_ntu has to converge to Q_con (m_ref*delta_h), we can safely - # calculate the output temperature. + T_in, T_sc, T_sh, T_out = get_condenser_phase_temperatures_and_alpha( + heat_exchanger=self, inputs=inputs, + Q_sc=Q_sc, Q_lat=Q_lat, Q_sh=Q_sh + ) - T_mean = inputs.T_con_in + self.calc_secondary_Q_flow(Q) / (self.m_flow_secondary_cp * 2) - tra_prop_med = self.calc_transport_properties_secondary_medium(T_mean) + tra_prop_med = self.calc_transport_properties_secondary_medium((T_in + T_out) / 2) alpha_med_wall = self.calc_alpha_secondary(tra_prop_med) - # Calculate secondary_medium side temperatures: - # Assumption loss is the same correlation for each regime - T_sc = inputs.T_con_in + self.calc_secondary_Q_flow(Q_sc) / self.m_flow_secondary_cp - T_sh = T_sc + self.calc_secondary_Q_flow(Q_lat) / self.m_flow_secondary_cp - T_out = T_sh + self.calc_secondary_Q_flow(Q_sh) / self.m_flow_secondary_cp - # 1. Regime: Subcooling Q_sc_ntu, A_sc = 0, 0 if Q_sc > 0 and (state_q0.T != self.state_outlet.T): @@ -165,13 +122,14 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_con) # Only use still available area: - A_sc = self.iterate_area(dT_max=(state_q0.T - inputs.T_con_in), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_sc) + A_sc = iterate_area(heat_exchanger=self, + dT_max=(state_q0.T - T_in), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_sc) A_sc = min(self.A, A_sc) - Q_sc_ntu, k_sc = self.calc_Q_ntu(dT_max=(state_q0.T - inputs.T_con_in), + Q_sc_ntu, k_sc = self.calc_Q_ntu(dT_max=(state_q0.T - T_in), alpha_pri=alpha_ref_wall, alpha_sec=alpha_med_wall, A=A_sc) @@ -188,10 +146,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): inputs=inputs ) - A_lat = self.iterate_area(dT_max=(state_q1.T - T_sc), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_lat) + A_lat = iterate_area(heat_exchanger=self, + dT_max=(state_q1.T - T_sc), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_lat) # Only use still available area: A_lat = min(self.A - A_sc, A_lat) @@ -221,20 +180,22 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): Q_ntu = Q_sh_ntu + Q_sc_ntu + Q_lat_ntu error = (Q_ntu / Q - 1) * 100 # Get possible dT_min: - dT_min_in = self.state_outlet.T - inputs.T_con_in + dT_min_in = self.state_outlet.T - T_in dT_min_out = self.state_inlet.T - T_out dT_min_LatSH = state_q1.T - T_sh - fs_state.set(name="A_con_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in condenser") + fs_state.set(name="A_con_sh", value=A_sh, unit="m2", + description="Area for superheat heat exchange in condenser") fs_state.set(name="A_con_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in condenser") - fs_state.set(name="A_con_sc", value=A_sc, unit="m2", description="Area for subcooling heat exchange in condenser") + fs_state.set(name="A_con_sc", value=A_sc, unit="m2", + description="Area for subcooling heat exchange in condenser") return error, min(dT_min_in, dT_min_LatSH, dT_min_out) -class MovingBoundaryNTUEvaporator(MovingBoundaryNTU): +class MovingBoundaryNTUEvaporator(BasicNTU): """ Evaporator class which implements the actual `calc` method. @@ -267,10 +228,12 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): self.calc_secondary_cp(T=inputs.T_eva_in) # First we separate the flow: - Q_sc, Q_lat, Q_sh, state_q0, state_q1 = self.separate_phases( - self.state_outlet, - self.state_inlet, - self.state_inlet.p + Q_sc, Q_lat, Q_sh, state_q0, state_q1 = separate_phases( + m_flow=self.m_flow, + med_prop=self.med_prop, + state_max=self.state_outlet, + state_min=self.state_inlet, + p=self.state_outlet.p ) Q = Q_sc + Q_lat + Q_sh @@ -295,10 +258,11 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_eva) if Q_lat > 0: - A_sh = self.iterate_area(dT_max=(inputs.T_eva_in - state_q1.T), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_sh) + A_sh = iterate_area(heat_exchanger=self, + dT_max=(inputs.T_eva_in - state_q1.T), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_sh) else: # if only sh is present --> full area: A_sh = self.A @@ -326,10 +290,10 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): ) if Q_sc > 0: - A_lat = self.iterate_area(dT_max=(T_sh - self.state_inlet.T), - alpha_pri=alpha_ref_wall, - alpha_sec=alpha_med_wall, - Q=Q_lat) + A_lat = iterate_area(heat_exchanger=self, dT_max=(T_sh - self.state_inlet.T), + alpha_pri=alpha_ref_wall, + alpha_sec=alpha_med_wall, + Q=Q_lat) else: A_lat = self.A - A_sh @@ -363,7 +327,9 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): dT_min_in = inputs.T_eva_in - self.state_outlet.T dT_min_out = T_out - self.state_inlet.T - fs_state.set(name="A_eva_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in evaporator") - fs_state.set(name="A_eva_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in evaporator") + fs_state.set(name="A_eva_sh", value=A_sh, unit="m2", + description="Area for superheat heat exchange in evaporator") + fs_state.set(name="A_eva_lat", value=A_lat, unit="m2", + description="Area for latent heat exchange in evaporator") return error, min(dT_min_out, dT_min_in) diff --git a/vclibpy/components/heat_exchangers/utils.py b/vclibpy/components/heat_exchangers/utils.py new file mode 100644 index 0000000..eb003f0 --- /dev/null +++ b/vclibpy/components/heat_exchangers/utils.py @@ -0,0 +1,93 @@ +import numpy as np + +from vclibpy import Inputs +from vclibpy.components.heat_exchangers import HeatExchanger + +from vclibpy.media import ThermodynamicState + + +def calc_area(Q, k, Tm): + if Tm * k == 0: + return 1e10 # A really large number, but not np.inf to still be able to calculate with it + else: + return Q / (k * Tm) + + +def calc_mean_temperature( + T_hot_in: float, T_hot_out: float, + T_cold_in: float, T_cold_out: float, + flow_type: str = "counter" +): + if flow_type == "counter": + dT_A = T_hot_in - T_cold_out + dT_B = T_hot_out - T_cold_in + elif flow_type == "parallel": + dT_A = T_hot_in - T_cold_in + dT_B = T_hot_out - T_cold_out + else: + raise TypeError("Given flow_type is not supported yet") + if dT_B < 0 or dT_A < 0: + return 0 # Heat can't be transferred + if np.isclose(dT_B, 0) or np.isclose(dT_A, 0) or np.isclose(dT_B, dT_A): + return abs((dT_A + dT_B) / 2) + return (dT_A - dT_B) / np.log(dT_A / dT_B) + + +def separate_phases(m_flow, med_prop, state_max: ThermodynamicState, state_min: ThermodynamicState, p: float): + """ + Separates a flow with possible phase changes into three parts: + subcooling (sc), latent phase change (lat), and superheating (sh) + at the given pressure. + + Args: + state_max (ThermodynamicState): State with higher enthalpy. + state_min (ThermodynamicState): State with lower enthalpy. + p (float): Pressure of phase change. + + Returns: + Tuple[float, float, float, ThermodynamicState, ThermodynamicState]: + Q_sc: Heat for subcooling. + Q_lat: Heat for latent phase change. + Q_sh: Heat for superheating. + state_q0: State at vapor quality 0 and the given pressure. + state_q1: State at vapor quality 1 and the given pressure. + """ + # Get relevant states: + state_q0 = med_prop.calc_state("PQ", p, 0) + state_q1 = med_prop.calc_state("PQ", p, 1) + Q_sc = max(0.0, + min((state_q0.h - state_min.h), + (state_max.h - state_min.h))) * m_flow + Q_lat = max(0.0, + (min(state_max.h, state_q1.h) - + max(state_min.h, state_q0.h))) * m_flow + Q_sh = max(0.0, + min((state_max.h - state_q1.h), + (state_max.h - state_min.h))) * m_flow + return Q_sc, Q_lat, Q_sh, state_q0, state_q1 + + +def get_condenser_phase_temperatures_and_alpha( + inputs: Inputs, + Q_sc: float, Q_lat: float, Q_sh: float, + heat_exchanger: HeatExchanger +): + cp = heat_exchanger.m_flow_secondary_cp + Q = Q_sc + Q_lat + Q_sh + + # Calculate secondary_medium side temperatures: + # Assumption loss is the same correlation for each regime + + # Calculate secondary_medium side temperatures: + # Assumption loss is the same correlation for each regime + if inputs.uses_condenser_inlet: + T_in = inputs.T_con_in + T_sc = inputs.T_con_in + heat_exchanger.calc_secondary_Q_flow(Q_sc) / cp + T_sh = T_sc + heat_exchanger.calc_secondary_Q_flow(Q_lat) / cp + T_out = T_sh + heat_exchanger.calc_secondary_Q_flow(Q_sh) / cp + else: + T_out = inputs.T_con_out + T_sh = T_out - heat_exchanger.calc_secondary_Q_flow(Q_sh) / cp + T_sc = T_sh - heat_exchanger.calc_secondary_Q_flow(Q_lat) / cp + T_in = T_sc - heat_exchanger.calc_secondary_Q_flow(Q_sc) / cp + return T_in, T_sc, T_sh, T_out diff --git a/vclibpy/datamodels.py b/vclibpy/datamodels.py index c2da3fe..d734e7e 100644 --- a/vclibpy/datamodels.py +++ b/vclibpy/datamodels.py @@ -134,7 +134,7 @@ def convert_to_str_value_format(self, with_unit_and_description: bool) -> Dict[s def get_name(self): return ";".join([k + "=" + str(round(v.value, 3)).replace(".", "_") - for k, v in self.get_variables().items()]) + for k, v in self.get_variables().items() if v.value is not None]) class FlowsheetState(VariableContainer): @@ -172,6 +172,7 @@ def __init__( n: float = None, T_eva_in: float = None, T_con_in: float = None, + T_con_out: float = None, m_flow_eva: float = None, m_flow_con: float = None, dT_eva_superheating: float = None, @@ -182,10 +183,17 @@ def __init__( Initializes an Inputs object with parameters representing external conditions for the vapor compression cycle. + Note that some inputs presuppose each other. For example, you can not + specify both `T_con_in`, `T_con_out`, and `m_flow_con`, + as one needs to be free for the simulation. + The current algorithms require `m_flow_con` + and either `T_con_in` or `T_con_out`. + Args: n (float): Relative compressor speed between 0 and 1 (unit: -). T_eva_in (float): Secondary side evaporator inlet temperature (unit: K). T_con_in (float): Secondary side condenser inlet temperature (unit: K). + T_con_out (float): Secondary side condenser inlet temperature (unit: K). m_flow_eva (float): Secondary side evaporator mass flow rate (unit: kg/s). m_flow_con (float): Secondary side condenser mass flow rate (unit: kg/s). dT_eva_superheating (float): Super-heating after evaporator (unit: K). @@ -193,6 +201,13 @@ def __init__( T_ambient (float): Ambient temperature of the machine (unit: K). """ super().__init__() + if T_con_out is not None and T_con_in is not None: + raise ValueError("You can either specify condenser inlet or " + "outlet temperature, not both.") + if T_con_out is None and T_con_in is None: + raise ValueError("You must either specify condenser inlet or " + "outlet temperature.") + self.set( name="n", value=n, @@ -211,6 +226,12 @@ def __init__( unit="K", description="Secondary side condenser inlet temperature" ) + self.set( + name="T_con_out", + value=T_con_out, + unit="K", + description="Secondary side condenser outlet temperature" + ) self.set( name="m_flow_con", value=m_flow_con, @@ -243,3 +264,19 @@ def __init__( unit="K", description="Ambient temperature of machine" ) + + @property + def T_con(self): + """Returns either T_con_in or T_con_out, depending on what is specifiec.""" + if self.T_con_in is not None: + return self.T_con_in + return self.T_con_out + + @property + def uses_condenser_inlet(self): + """Indicates if condenser inlet or outlet is specified""" + return self.T_con_in is not None + + +if __name__ == '__main__': + print(Inputs().T_con_out) diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index 79621ea..cb046d7 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -73,13 +73,20 @@ def get_all_components(self) -> List[BaseComponent]: return [self.condenser, self.evaporator] def get_start_condensing_pressure(self, inputs: Inputs): - T_3_start = inputs.T_con_in + inputs.dT_con_subcooling + if inputs.uses_condenser_inlet: + T_3_start = inputs.T_con_in + inputs.dT_con_subcooling + else: + try: + dT_con_start = inputs.dT_con_start + except AttributeError: + dT_con_start = 10 + T_3_start = inputs.T_con_out - dT_con_start p_2_start = self.med_prop.calc_state("TQ", T_3_start, 0).p return p_2_start def improve_first_condensing_guess(self, inputs: Inputs, m_flow_guess, p_2_guess, dT_pinch_assumption=0): self.condenser.m_flow_secondary = inputs.m_flow_con # [kg/s] - self.condenser.calc_secondary_cp(T=inputs.T_con_in) + self.condenser.calc_secondary_cp(T=inputs.T_con) def nonlinear_func(p_2, *args): _flowsheet, _inputs, _m_flow_ref, _dT_pinch = args @@ -421,7 +428,12 @@ def calculate_outputs_for_valid_pressures( error_con, dT_min_con = self.condenser.calc(inputs=inputs, fs_state=fs_state) error_eva, dT_min_eva = self.evaporator.calc(inputs=inputs, fs_state=fs_state) P_el = self.calc_electrical_power(fs_state=fs_state, inputs=inputs) - T_con_out = inputs.T_con_in + Q_con_outer / self.condenser.m_flow_secondary_cp + if inputs.uses_condenser_inlet: + T_con_in = inputs.T_con_in + T_con_out = T_con_in + Q_con_outer / self.condenser.m_flow_secondary_cp + else: + T_con_out = inputs.T_con_out + T_con_in = T_con_out - Q_con_outer / self.condenser.m_flow_secondary_cp # COP based on P_el and Q_con: COP_inner = Q_con / P_el @@ -430,6 +442,17 @@ def calculate_outputs_for_valid_pressures( COP_carnot = (T_con_out / (T_con_out - inputs.T_eva_in)) carnot_quality = COP_inner / COP_carnot # Calc return temperature: + if inputs.uses_condenser_inlet: + fs_state.set( + name="T_con_out", value=T_con_out, unit="K", + description="Condenser outlet temperature" + ) + else: + fs_state.set( + name="T_con_in", value=T_con_in, unit="K", + description="Condenser inlet temperature" + ) + fs_state.set( name="P_el", value=P_el, unit="W", description="Power consumption" @@ -576,13 +599,19 @@ def _plot_secondary_heat_flow_rates(self, ax, inputs): self.evaporator.state_outlet.h * self.evaporator.m_flow - Q_eva ]) / self.evaporator.m_flow self.condenser.m_flow_secondary = inputs.m_flow_con - self.condenser.calc_secondary_cp(T=inputs.T_con_in) + self.condenser.calc_secondary_cp(T=inputs.T_con) self.evaporator.m_flow_secondary = inputs.m_flow_eva self.evaporator.calc_secondary_cp(T=inputs.T_eva_in) - ax.plot(delta_H_con / 1000, [ - inputs.T_con_in - 273.15, - inputs.T_con_in + Q_con / self.condenser.m_flow_secondary_cp - 273.15 - ], color="b") + if inputs.uses_condenser_inlet: + ax.plot(delta_H_con / 1000, [ + inputs.T_con_in - 273.15, + inputs.T_con_in + Q_con / self.condenser.m_flow_secondary_cp - 273.15 + ], color="b") + else: + ax.plot(delta_H_con / 1000, [ + inputs.T_con_out - Q_con / self.condenser.m_flow_secondary_cp - 273.15, + inputs.T_con_out - 273.15 + ], color="b") ax.plot(delta_H_eva / 1000, [ inputs.T_eva_in - 273.15, inputs.T_eva_in - Q_eva / self.evaporator.m_flow_secondary_cp - 273.15 diff --git a/vclibpy/utils/automation.py b/vclibpy/utils/automation.py index 2885b01..2578d2a 100644 --- a/vclibpy/utils/automation.py +++ b/vclibpy/utils/automation.py @@ -55,15 +55,18 @@ def calc_multiple_states( def full_factorial_map_generation( heat_pump: BaseCycle, T_eva_in_ar: Union[list, np.ndarray], - T_con_in_ar: Union[list, np.ndarray], + T_con_ar: Union[list, np.ndarray], n_ar: Union[list, np.ndarray], m_flow_con: float, m_flow_eva: float, save_path: Union[pathlib.Path, str], dT_eva_superheating=5, dT_con_subcooling=0, + use_condenser_inlet: bool = True, + dT_con_start: float = 10, use_multiprocessing: bool = False, save_plots: bool = False, + raise_errors: bool = False, **kwargs ) -> (pathlib.Path, pathlib.Path): """ @@ -71,24 +74,39 @@ def full_factorial_map_generation( used in other simulation tools like Modelica or to analyze the off-design of the flowsheet. The results are stored and returned as .sdf and .csv files. - Currently, only varying T_eva_in, T_con_in, and n is implemented. + Currently, only varying T_eva_in, T_con_in (or T_con_out), and n is implemented. However, changing this to more dimensions or other variables is not much work. In this case, please raise an issue. Args: heat_pump (BaseCycle): The flowsheet to use - T_eva_in_ar: Array with inputs for T_eva_in - T_con_in_ar: Array with inputs for T_con_in - n_ar: Array with inputs for n_ar - m_flow_con: Condenser mass flow rate - m_flow_eva: Evaporator mass flow rate - save_path: Where to save all results. - dT_eva_superheating: Evaporator superheating - dT_con_subcooling: Condenser subcooling + T_eva_in_ar (list): + Array with inputs for T_eva_in + T_con_ar (list): + Array with inputs for T_con_in or T_con_out, see `use_condenser_inlet` + n_ar (list): + Array with inputs for n_ar + m_flow_con (float): + Condenser mass flow rate + m_flow_eva (float): + Evaporator mass flow rate + save_path (Path): + Where to save all results. + dT_eva_superheating (float): + Evaporator superheating + dT_con_subcooling (float): + Condenser subcooling + use_condenser_inlet (bool): + True to consider T_con_ar as inlet, false for outlet. + dT_con_start (float): + Guess value for starting condenser pressure, + only relevant if `use_condenser_inlet=False` use_multiprocessing: True to use multiprocessing. May speed up the calculation. Default is False save_plots: True to save plots of each steady state point. Default is False + raise_errors: + True to raise errors if they occur. **kwargs: Solver settings for the flowsheet Returns: @@ -106,16 +124,28 @@ def full_factorial_map_generation( idx_for_access_later = [] for i_T_eva_in, T_eva_in in enumerate(T_eva_in_ar): for i_n, n in enumerate(n_ar): - for i_T_con_in, T_con_in in enumerate(T_con_in_ar): - idx_for_access_later.append([i_n, i_T_con_in, i_T_eva_in]) - inputs = Inputs(n=n, - T_eva_in=T_eva_in, - T_con_in=T_con_in, - m_flow_eva=m_flow_eva, - m_flow_con=m_flow_con, - dT_eva_superheating=dT_eva_superheating, - dT_con_subcooling=dT_con_subcooling) - list_mp_inputs.append([heat_pump, inputs, kwargs]) + for i_T_con, T_con in enumerate(T_con_ar): + idx_for_access_later.append([i_n, i_T_con, i_T_eva_in]) + base_inputs = dict( + n=n, + T_eva_in=T_eva_in, + m_flow_eva=m_flow_eva, + m_flow_con=m_flow_con, + dT_eva_superheating=dT_eva_superheating, + dT_con_subcooling=dT_con_subcooling + ) + if use_condenser_inlet: + inputs = Inputs( + T_con_in=T_con, + **base_inputs + ) + else: + inputs = Inputs( + T_con_out=T_con, + **base_inputs + ) + inputs.set(name="dT_con_start", value=dT_con_start, unit="K") + list_mp_inputs.append([heat_pump, inputs, kwargs, raise_errors]) list_inputs.append(inputs) fs_states = [] i = 0 @@ -127,13 +157,13 @@ def full_factorial_map_generation( logger.info(f"Ran {i} of {len(list_mp_inputs)} points") else: for inputs in list_inputs: - fs_state = _calc_single_hp_state([heat_pump, inputs, kwargs]) + fs_state = _calc_single_hp_state([heat_pump, inputs, kwargs, raise_errors]) fs_states.append(fs_state) i += 1 logger.info(f"Ran {i} of {len(list_mp_inputs)} points") # Save to sdf - result_shape = (len(n_ar), len(T_con_in_ar), len(T_eva_in_ar)) + result_shape = (len(n_ar), len(T_con_ar), len(T_eva_in_ar)) _dummy = np.zeros(result_shape) # Use a copy to avoid overwriting of values of any sort. _dummy[:] = np.nan # Get all possible values: @@ -156,9 +186,9 @@ def full_factorial_map_generation( ) for fs_state, idx_triple in zip(fs_states, idx_for_access_later): - i_n, i_T_con_in, i_T_eva_in = idx_triple + i_n, i_T_con, i_T_eva_in = idx_triple for variable_name, variable in fs_state.get_variables().items(): - all_variables[variable_name][i_n][i_T_con_in][i_T_eva_in] = variable.value + all_variables[variable_name][i_n][i_T_con][i_T_eva_in] = variable.value _nd_data = {} for variable, nd_data in all_variables.items(): @@ -171,7 +201,7 @@ def full_factorial_map_generation( _scale_values = { "n": n_ar, - "T_con_in": T_con_in_ar, + "T_con_in" if use_condenser_inlet else "T_con_out": T_con_ar, "T_eva_in": T_eva_in_ar } inputs: Inputs = list_inputs[0] @@ -207,12 +237,14 @@ def full_factorial_map_generation( def _calc_single_hp_state(data): """Helper function for a single state to enable multiprocessing""" - heat_pump, inputs, kwargs = data + heat_pump, inputs, kwargs, raise_errors = data fs_state = None try: fs_state = heat_pump.calc_steady_state(inputs=inputs, **kwargs) except Exception as e: + if raise_errors: + raise e logger.error(f"An error occurred for input: {inputs.__dict__}: {e}") if fs_state is None: fs_state = FlowsheetState() From 81a0ad2caccbebe7b4abf2b774403d8a002b3ff1 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Wed, 15 May 2024 12:21:27 +0200 Subject: [PATCH 11/36] fix examples and add case for T_con_out in examples --- examples/e6_simple_heat_pump.py | 11 ++++++----- tests/test_examples.py | 7 +++++-- vclibpy/datamodels.py | 3 --- vclibpy/utils/automation.py | 2 +- 4 files changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/e6_simple_heat_pump.py b/examples/e6_simple_heat_pump.py index 7714f7e..9354d9c 100644 --- a/examples/e6_simple_heat_pump.py +++ b/examples/e6_simple_heat_pump.py @@ -1,6 +1,6 @@ # # Example for a heat pump with a standard cycle -def main(): +def main(use_condenser_inlet: bool = True): # Let's start the complete cycle simulation with the # most basic flowsheet, the standard-cycle. As all flowsheets # contain a condenser and an evaporator, we defined a common BaseCycle @@ -62,7 +62,7 @@ def main(): # As in the other example, we can specify save-paths, # solver settings and inputs to vary: # Note that T_con can either be inlet or outlet, depending on the setting - # of `use_condenser_inlet` later on. In this case, we simulate the inlet, T_con_in + # of `use_condenser_inlet`. Per default, we simulate the inlet, T_con_in save_path = r"D:\00_temp\simple_heat_pump" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] T_con_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] @@ -82,7 +82,7 @@ def main(): T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, - use_condenser_inlet=True, + use_condenser_inlet=use_condenser_inlet, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -113,16 +113,17 @@ def main(): # We can also use existing 3D-plotting scripts in vclibpy to analyze the # dependencies. For this, we need the .sdf file. In the sdf, the field names are without # the unit and description, as those are accessible in the file-format in other columns. + # Depending on whether we varied the inlet or outlet, we have to specify the correct name. from vclibpy.utils.plotting import plot_sdf_map plot_sdf_map( filepath_sdf=save_path_sdf, nd_data="COP", first_dimension="T_eva_in", - second_dimension="T_con_in", + second_dimension="T_con_in" if use_condenser_inlet else "T_con_out", fluids=["Propane"], flowsheets=["Standard"] ) if __name__ == "__main__": - main() + main(use_condenser_inlet=False) diff --git a/tests/test_examples.py b/tests/test_examples.py index cb8953d..1dde3ad 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -13,12 +13,12 @@ class TestExamples(unittest.TestCase): def setUp(self) -> None: self.timeout = 100 # Seconds which the script is allowed to run - def _run_example(self, example, timeout=None): + def _run_example(self, example, timeout=None, **kwargs): ex_py = pathlib.Path(__file__).absolute().parents[1].joinpath("examples", example) sys.path.insert(0, str(ex_py.parent)) module = importlib.import_module(ex_py.stem) example_main = getattr(module, "main") - example_main() + example_main(**kwargs) def test_e1_refrigerant_data(self): self._run_example(example="e1_refrigerant_data.py") @@ -38,6 +38,9 @@ def test_e5_expansion_valve(self): def test_e6_simple_heat_pump(self): self._run_example(example="e6_simple_heat_pump.py") + def test_e6_simple_heat_pump_outlet(self): + self._run_example(example="e6_simple_heat_pump.py", use_condenser_inlet=False) + def test_e7_vapor_injection(self): self._run_example(example="e7_vapor_injection.py") diff --git a/vclibpy/datamodels.py b/vclibpy/datamodels.py index d734e7e..6a0c5a4 100644 --- a/vclibpy/datamodels.py +++ b/vclibpy/datamodels.py @@ -204,9 +204,6 @@ def __init__( if T_con_out is not None and T_con_in is not None: raise ValueError("You can either specify condenser inlet or " "outlet temperature, not both.") - if T_con_out is None and T_con_in is None: - raise ValueError("You must either specify condenser inlet or " - "outlet temperature.") self.set( name="n", diff --git a/vclibpy/utils/automation.py b/vclibpy/utils/automation.py index 2578d2a..3826ca8 100644 --- a/vclibpy/utils/automation.py +++ b/vclibpy/utils/automation.py @@ -207,7 +207,7 @@ def full_factorial_map_generation( inputs: Inputs = list_inputs[0] _parameters = {} for name, variable in inputs.items(): - if name not in _scale_values: + if name not in list(_scale_values.keys()) + ["T_con_out", "T_con_in"]: _parameters[name] = { "data": variable.value, "unit": variable.unit, From a921fdff9bd21733635b44c919e16e15872c5331 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Wed, 15 May 2024 13:57:21 +0200 Subject: [PATCH 12/36] handle None inputs --- vclibpy/datamodels.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/vclibpy/datamodels.py b/vclibpy/datamodels.py index 6a0c5a4..a77ae7d 100644 --- a/vclibpy/datamodels.py +++ b/vclibpy/datamodels.py @@ -129,12 +129,17 @@ def convert_to_str_value_format(self, with_unit_and_description: bool) -> Dict[s dict: Containing all variables and values. """ if with_unit_and_description: - return {f"{k} in {v.unit} ({v.description})": v.value for k, v in self.items()} - return {k: v.value for k, v in self.items()} + return {f"{k} in {v.unit} ({v.description})": v.value for k, v in self.items_not_none()} + return {k: v.value for k, v in self.items_not_none()} def get_name(self): + """Get the name based on variable names and rounded values""" return ";".join([k + "=" + str(round(v.value, 3)).replace(".", "_") - for k, v in self.get_variables().items() if v.value is not None]) + for k, v in self.items_not_none()]) + + def items_not_none(self): + """Get all items where the value is not None""" + return {k: v for k, v in self.items() if v.value is not None}.items() class FlowsheetState(VariableContainer): From 4c904807d5935f00e23277307be4130ae4e0e0a7 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Wed, 15 May 2024 13:57:36 +0200 Subject: [PATCH 13/36] avoid runtime error for floating precision errors --- .../heat_exchangers/moving_boundary_Tm.py | 14 ++++++++------ .../heat_exchangers/moving_boundary_ntu.py | 12 ++++++------ 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py index 649d432..69c301e 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py @@ -1,5 +1,7 @@ import logging +import numpy as np + from vclibpy.components.heat_exchangers import utils from vclibpy.datamodels import FlowsheetState, Inputs from vclibpy.components.heat_exchangers import HeatExchanger @@ -59,7 +61,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 1. Regime: Subcooling Q_sc_Tm, A_sc, A_sc_available = 0, 0, 0 - if Q_sc > 0 and (state_q0.T != self.state_outlet.T): + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_outlet.T): # Get transport properties: tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(state_q0, self.state_outlet) alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_con) @@ -76,7 +78,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 2. Regime: Latent heat exchange Q_lat_Tm, A_lat, A_lat_available = 0, 0, 0 - if Q_lat > 0: + if not np.isclose(Q_lat, 0): # Get transport properties: alpha_ref_wall = self.calc_alpha_two_phase( state_q0=state_q0, @@ -98,7 +100,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 3. Regime: Superheat heat exchange Q_sh_Tm, A_sh = 0, 0 - if Q_sh and (self.state_inlet.T != state_q1.T): + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_inlet.T, state_q1.T): # Get transport properties: tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(self.state_inlet, state_q1) alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_con) @@ -190,7 +192,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 1. Regime: Superheating Q_sh_Tm, A_sh, A_sh_available = 0, 0, 0 - if Q_sh and (self.state_outlet.T != state_q1.T): + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_outlet.T, state_q1.T): # Get transport properties: tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(self.state_outlet, state_q1) alpha_ref_wall = self.calc_alpha_gas(tra_prop_ref_eva) @@ -209,7 +211,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 2. Regime: Latent heat exchange Q_lat_Tm, A_lat, A_lat_available = 0, 0, 0 - if Q_lat > 0: + if not np.isclose(Q_lat, 0): alpha_ref_wall = self.calc_alpha_two_phase( state_q0=state_q0, state_q1=state_q1, @@ -230,7 +232,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 3. Regime: Subcooling Q_sc_Tm, A_sc = 0, 0 - if Q_sc > 0 and (state_q0.T != self.state_inlet.T): + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_inlet.T): # Get transport properties: tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(state_q0, self.state_inlet) alpha_ref_wall = self.calc_alpha_liquid(tra_prop_ref_eva) diff --git a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py index 4f68d8b..e1cbe14 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_ntu.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_ntu.py @@ -115,7 +115,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 1. Regime: Subcooling Q_sc_ntu, A_sc = 0, 0 - if Q_sc > 0 and (state_q0.T != self.state_outlet.T): + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_outlet.T): self.set_primary_cp((state_q0.h - self.state_outlet.h) / (state_q0.T - self.state_outlet.T)) # Get transport properties: tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(state_q0, self.state_outlet) @@ -136,7 +136,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 2. Regime: Latent heat exchange Q_lat_ntu, A_lat = 0, 0 - if Q_lat > 0: + if not np.isclose(Q_lat, 0): self.set_primary_cp(np.inf) # Get transport properties: alpha_ref_wall = self.calc_alpha_two_phase( @@ -162,7 +162,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 3. Regime: Superheat heat exchange Q_sh_ntu, A_sh = 0, 0 - if Q_sh and (self.state_inlet.T != state_q1.T): + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_inlet.T, state_q1.T): self.set_primary_cp((self.state_inlet.h - state_q1.h) / (self.state_inlet.T - state_q1.T)) # Get transport properties: tra_prop_ref_con = self.med_prop.calc_mean_transport_properties(self.state_inlet, state_q1) @@ -251,7 +251,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 1. Regime: Superheating Q_sh_ntu, A_sh = 0, 0 - if Q_sh and (self.state_outlet.T != state_q1.T): + if not np.isclose(Q_sh, 0) and not np.isclose(self.state_outlet.T, state_q1.T): self.set_primary_cp((self.state_outlet.h - state_q1.h) / (self.state_outlet.T - state_q1.T)) # Get transport properties: tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(self.state_outlet, state_q1) @@ -279,7 +279,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 2. Regime: Latent heat exchange Q_lat_ntu, A_lat = 0, 0 - if Q_lat > 0: + if not np.isclose(Q_lat, 0): self.set_primary_cp(np.inf) alpha_ref_wall = self.calc_alpha_two_phase( @@ -307,7 +307,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): # 3. Regime: Subcooling Q_sc_ntu, A_sc = 0, 0 - if Q_sc > 0 and (state_q0.T != self.state_inlet.T): + if not np.isclose(Q_sc, 0) and not np.isclose(state_q0.T, self.state_inlet.T): self.set_primary_cp((state_q0.h - self.state_inlet.h) / (state_q0.T - self.state_inlet.T)) # Get transport properties: tra_prop_ref_eva = self.med_prop.calc_mean_transport_properties(state_q0, self.state_inlet) From a126fc0738320af016d6047f85d59e8ae2f45b50 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Wed, 15 May 2024 16:32:09 +0200 Subject: [PATCH 14/36] print for CI debugging --- examples/e6_simple_heat_pump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/e6_simple_heat_pump.py b/examples/e6_simple_heat_pump.py index 9354d9c..ce9b1f9 100644 --- a/examples/e6_simple_heat_pump.py +++ b/examples/e6_simple_heat_pump.py @@ -98,7 +98,7 @@ def main(use_condenser_inlet: bool = True): # the values using e.g. pandas. It is also the second return value of the function. import pandas as pd df = pd.read_csv(save_path_csv, index_col=0) - df + print(df) # Now, we can plot variables, for example as a scatter plot using matplotlib. # You have to know the names, which are the column headers. import matplotlib.pyplot as plt From 727137663daa1a8162e39d2981622d6bd101668d Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Wed, 15 May 2024 17:38:08 +0200 Subject: [PATCH 15/36] specify s for debugging --- examples/e6_simple_heat_pump.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/e6_simple_heat_pump.py b/examples/e6_simple_heat_pump.py index ce9b1f9..48854a6 100644 --- a/examples/e6_simple_heat_pump.py +++ b/examples/e6_simple_heat_pump.py @@ -98,13 +98,13 @@ def main(use_condenser_inlet: bool = True): # the values using e.g. pandas. It is also the second return value of the function. import pandas as pd df = pd.read_csv(save_path_csv, index_col=0) - print(df) + df # Now, we can plot variables, for example as a scatter plot using matplotlib. # You have to know the names, which are the column headers. import matplotlib.pyplot as plt x_name = "n in - (Relative compressor speed)" y_name = "COP in - (Coefficient of performance)" - plt.scatter(df[x_name], df[y_name]) + plt.scatter(df[x_name], df[y_name], s=20) plt.ylabel(y_name) plt.xlabel(x_name) plt.show() From 53e86921dacec57154e4a0630e2df12c1885ecf7 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst@eonerc.rwth-aachen.de" Date: Wed, 15 May 2024 15:47:21 +0000 Subject: [PATCH 16/36] chore(examples): Automatic commit of example files in Markdown and Jupyter Notebook format. --- docs/jupyter_notebooks/e6_simple_heat_pump.ipynb | 12 ++++++------ docs/jupyter_notebooks/e7_vapor_injection.ipynb | 6 +++--- docs/source/examples/e6_simple_heat_pump.md | 12 ++++++++---- docs/source/examples/e7_vapor_injection.md | 7 ++++--- 4 files changed, 21 insertions(+), 16 deletions(-) diff --git a/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb b/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb index ea8dfc8..da02875 100644 --- a/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb +++ b/docs/jupyter_notebooks/e6_simple_heat_pump.ipynb @@ -44,14 +44,14 @@ { "cell_type": "markdown", "metadata": {}, - "source": "As in the other example, we can specify save-paths,\nsolver settings and inputs to vary:\n" + "source": "As in the other example, we can specify save-paths,\nsolver settings and inputs to vary:\nNote that T_con can either be inlet or outlet, depending on the setting\nof `use_condenser_inlet`. Per default, we simulate the inlet, T_con_in\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": "save_path = r\"D:\\00_temp\\simple_heat_pump\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_in_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" + "source": "save_path = r\"D:\\00_temp\\simple_heat_pump\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" }, { "cell_type": "markdown", @@ -63,7 +63,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nsave_path_sdf, save_path_csv = utils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_in_ar=T_con_in_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" + "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nsave_path_sdf, save_path_csv = utils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_ar=T_con_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_condenser_inlet=use_condenser_inlet,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" }, { "cell_type": "markdown", @@ -99,19 +99,19 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "import matplotlib.pyplot as plt\nx_name = \"n in - (Relative compressor speed)\"\ny_name = \"COP in - (Coefficient of performance)\"\nplt.scatter(df[x_name], df[y_name])\nplt.ylabel(y_name)\nplt.xlabel(x_name)\nplt.show()\n" + "source": "import matplotlib.pyplot as plt\nx_name = \"n in - (Relative compressor speed)\"\ny_name = \"COP in - (Coefficient of performance)\"\nplt.scatter(df[x_name], df[y_name], s=20)\nplt.ylabel(y_name)\nplt.xlabel(x_name)\nplt.show()\n" }, { "cell_type": "markdown", "metadata": {}, - "source": "Looking at the results, we see that a higher frequency often leads to lower COP values.\nHowever, other inputs (temperatures) have a greater impact on the COP.\nWe can also use existing 3D-plotting scripts in vclibpy to analyze the\ndependencies. For this, we need the .sdf file. In the sdf, the field names are without\nthe unit and description, as those are accessible in the file-format in other columns.\n" + "source": "Looking at the results, we see that a higher frequency often leads to lower COP values.\nHowever, other inputs (temperatures) have a greater impact on the COP.\nWe can also use existing 3D-plotting scripts in vclibpy to analyze the\ndependencies. For this, we need the .sdf file. In the sdf, the field names are without\nthe unit and description, as those are accessible in the file-format in other columns.\nDepending on whether we varied the inlet or outlet, we have to specify the correct name.\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": "from vclibpy.utils.plotting import plot_sdf_map\nplot_sdf_map(\n filepath_sdf=save_path_sdf,\n nd_data=\"COP\",\n first_dimension=\"T_eva_in\",\n second_dimension=\"T_con_in\",\n fluids=[\"Propane\"],\n flowsheets=[\"Standard\"]\n)\n" + "source": "from vclibpy.utils.plotting import plot_sdf_map\nplot_sdf_map(\n filepath_sdf=save_path_sdf,\n nd_data=\"COP\",\n first_dimension=\"T_eva_in\",\n second_dimension=\"T_con_in\" if use_condenser_inlet else \"T_con_out\",\n fluids=[\"Propane\"],\n flowsheets=[\"Standard\"]\n)\n" } ], "metadata": { diff --git a/docs/jupyter_notebooks/e7_vapor_injection.ipynb b/docs/jupyter_notebooks/e7_vapor_injection.ipynb index f2599e1..2519f73 100644 --- a/docs/jupyter_notebooks/e7_vapor_injection.ipynb +++ b/docs/jupyter_notebooks/e7_vapor_injection.ipynb @@ -63,7 +63,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "save_path = r\"D:\\00_temp\\vapor_injection\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_in_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" + "source": "save_path = r\"D:\\00_temp\\vapor_injection\"\nT_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15]\nT_con_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15]\nn_ar = [0.3, 0.7, 1]\n" }, { "cell_type": "markdown", @@ -75,7 +75,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_in_ar=T_con_in_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" + "source": "import logging\nlogging.basicConfig(level=\"INFO\")\n\nfrom vclibpy import utils\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=save_path,\n T_con_ar=T_con_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_condenser_inlet=True,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" }, { "cell_type": "markdown", @@ -123,7 +123,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "heat_pump = VaporInjectionEconomizer(\n evaporator=evaporator,\n condenser=condenser,\n fluid=\"Propane\",\n economizer=economizer,\n high_pressure_compressor=high_pressure_compressor,\n low_pressure_compressor=low_pressure_compressor,\n high_pressure_valve=high_pressure_valve,\n low_pressure_valve=low_pressure_valve\n)\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=r\"D:\\00_temp\\vapor_injection_economizer\",\n T_con_in_ar=T_con_in_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" + "source": "heat_pump = VaporInjectionEconomizer(\n evaporator=evaporator,\n condenser=condenser,\n fluid=\"Propane\",\n economizer=economizer,\n high_pressure_compressor=high_pressure_compressor,\n low_pressure_compressor=low_pressure_compressor,\n high_pressure_valve=high_pressure_valve,\n low_pressure_valve=low_pressure_valve\n)\nutils.full_factorial_map_generation(\n heat_pump=heat_pump,\n save_path=r\"D:\\00_temp\\vapor_injection_economizer\",\n T_con_ar=T_con_ar,\n T_eva_in_ar=T_eva_in_ar,\n n_ar=n_ar,\n use_multiprocessing=False,\n save_plots=True,\n m_flow_con=0.2,\n m_flow_eva=0.9,\n dT_eva_superheating=5,\n dT_con_subcooling=0,\n)\n" } ], "metadata": { diff --git a/docs/source/examples/e6_simple_heat_pump.md b/docs/source/examples/e6_simple_heat_pump.md index 7b5c84e..988de47 100644 --- a/docs/source/examples/e6_simple_heat_pump.md +++ b/docs/source/examples/e6_simple_heat_pump.md @@ -71,11 +71,13 @@ heat_pump = StandardCycle( As in the other example, we can specify save-paths, solver settings and inputs to vary: +Note that T_con can either be inlet or outlet, depending on the setting +of `use_condenser_inlet`. Per default, we simulate the inlet, T_con_in ```python save_path = r"D:\00_temp\simple_heat_pump" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] -T_con_in_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] +T_con_ar = [30 + 273.15, 50 + 273.15, 70 + 273.15] n_ar = [0.3, 0.7, 1] ``` @@ -92,9 +94,10 @@ from vclibpy import utils save_path_sdf, save_path_csv = utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=use_condenser_inlet, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -127,7 +130,7 @@ You have to know the names, which are the column headers. import matplotlib.pyplot as plt x_name = "n in - (Relative compressor speed)" y_name = "COP in - (Coefficient of performance)" -plt.scatter(df[x_name], df[y_name]) +plt.scatter(df[x_name], df[y_name], s=20) plt.ylabel(y_name) plt.xlabel(x_name) plt.show() @@ -138,6 +141,7 @@ However, other inputs (temperatures) have a greater impact on the COP. We can also use existing 3D-plotting scripts in vclibpy to analyze the dependencies. For this, we need the .sdf file. In the sdf, the field names are without the unit and description, as those are accessible in the file-format in other columns. +Depending on whether we varied the inlet or outlet, we have to specify the correct name. ```python from vclibpy.utils.plotting import plot_sdf_map @@ -145,7 +149,7 @@ plot_sdf_map( filepath_sdf=save_path_sdf, nd_data="COP", first_dimension="T_eva_in", - second_dimension="T_con_in", + second_dimension="T_con_in" if use_condenser_inlet else "T_con_out", fluids=["Propane"], flowsheets=["Standard"] ) diff --git a/docs/source/examples/e7_vapor_injection.md b/docs/source/examples/e7_vapor_injection.md index 05e47cf..26d96c1 100644 --- a/docs/source/examples/e7_vapor_injection.md +++ b/docs/source/examples/e7_vapor_injection.md @@ -85,7 +85,7 @@ solver settings and inputs to vary: ```python save_path = r"D:\00_temp\vapor_injection" T_eva_in_ar = [-10 + 273.15, 273.15, 10 + 273.15] -T_con_in_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] +T_con_ar = [30 + 273.15, 50 + 273.15, 60 + 273.15] n_ar = [0.3, 0.7, 1] ``` @@ -102,9 +102,10 @@ from vclibpy import utils utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=save_path, - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, + use_condenser_inlet=True, use_multiprocessing=False, save_plots=True, m_flow_con=0.2, @@ -161,7 +162,7 @@ heat_pump = VaporInjectionEconomizer( utils.full_factorial_map_generation( heat_pump=heat_pump, save_path=r"D:\00_temp\vapor_injection_economizer", - T_con_in_ar=T_con_in_ar, + T_con_ar=T_con_ar, T_eva_in_ar=T_eva_in_ar, n_ar=n_ar, use_multiprocessing=False, From 619c1e371ba0dec3696be77b275c8303535774df Mon Sep 17 00:00:00 2001 From: Julius Klupp Date: Sat, 25 May 2024 11:36:04 +0200 Subject: [PATCH 17/36] adjustments for multiprocessing for bes-rules compatability --- vclibpy/components/compressors/compressor.py | 4 ++-- vclibpy/flowsheets/base.py | 3 ++- vclibpy/utils/automation.py | 11 ++++++----- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/vclibpy/components/compressors/compressor.py b/vclibpy/components/compressors/compressor.py index 35d2a74..e4bb261 100644 --- a/vclibpy/components/compressors/compressor.py +++ b/vclibpy/components/compressors/compressor.py @@ -127,7 +127,7 @@ def calc_state_outlet(self, p_outlet: float, inputs: Inputs, fs_state: Flowsheet self.state_inlet.h + (state_outlet_isentropic.h - self.state_inlet.h) / eta_is ) - fs_state.set(name="eta_is", value=eta_is, unit="%", description="Isentropic efficiency") + fs_state.set(name="eta_is", value=eta_is, unit="-", description="Isentropic efficiency") self.state_outlet = self.med_prop.calc_state("PH", p_outlet, h_outlet) def calc_m_flow(self, inputs: Inputs, fs_state: FlowsheetState) -> float: @@ -148,7 +148,7 @@ def calc_m_flow(self, inputs: Inputs, fs_state: FlowsheetState) -> float: self.get_n_absolute(inputs.n) ) self.m_flow = self.state_inlet.d * V_flow_ref - fs_state.set(name="lambda_h", value=lambda_h, unit="%", description="Volumetric efficiency") + fs_state.set(name="lambda_h", value=lambda_h, unit="-", description="Volumetric efficiency") fs_state.set(name="V_flow_ref", value=V_flow_ref, unit="m3/s", description="Refrigerant volume flow rate") fs_state.set(name="m_flow_ref", value=self.m_flow, unit="kg/s", description="Refrigerant mass flow rate") return self.m_flow diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index cb046d7..6a4eda2 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -65,7 +65,8 @@ def setup_new_fluid(self, fluid): self.fluid = fluid def terminate(self): - self.med_prop.terminate() + if self.med_prop is not None: + self.med_prop.terminate() for component in self.get_all_components(): component.terminate_secondary_med_prop() diff --git a/vclibpy/utils/automation.py b/vclibpy/utils/automation.py index 3826ca8..02a0b85 100644 --- a/vclibpy/utils/automation.py +++ b/vclibpy/utils/automation.py @@ -150,7 +150,8 @@ def full_factorial_map_generation( fs_states = [] i = 0 if use_multiprocessing: - pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) + # pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) + pool = multiprocessing.Pool(processes=10) for fs_state in pool.imap(_calc_single_hp_state, list_mp_inputs): fs_states.append(fs_state) i += 1 @@ -174,15 +175,15 @@ def full_factorial_map_generation( all_variables.update({var: _dummy.copy() for var in fs_state.get_variable_names()}) all_variables_info.update({var: variable for var, variable in fs_state.get_variables().items()}) variables_to_excel.append({ - **inputs.convert_to_str_value_format(with_unit_and_description=True), - **fs_state.convert_to_str_value_format(with_unit_and_description=True), + **inputs.convert_to_str_value_format(with_unit_and_description=False), + **fs_state.convert_to_str_value_format(with_unit_and_description=False), }) # Save to excel save_path_sdf = save_path.joinpath(f"{heat_pump.flowsheet_name}_{heat_pump.fluid}.sdf") save_path_csv = save_path.joinpath(f"{heat_pump.flowsheet_name}_{heat_pump.fluid}.csv") pd.DataFrame(variables_to_excel).to_csv( - save_path_csv + save_path_csv, sep=";" ) for fs_state, idx_triple in zip(fs_states, idx_for_access_later): @@ -227,7 +228,7 @@ def full_factorial_map_generation( heat_pump.fluid: (_scales, _nd_data, _parameters) } } - utils.save_to_sdf(data=sdf_data, save_path=save_path_sdf) + # utils.save_to_sdf(data=sdf_data, save_path=save_path_sdf) # Terminate heat pump med-props: heat_pump.terminate() From 7bf99eeadbdf89c47df8b3b2ebd0623ddc08d2be Mon Sep 17 00:00:00 2001 From: Julius Klupp Date: Mon, 10 Jun 2024 16:16:13 +0200 Subject: [PATCH 18/36] Bug fixes on TenCoefficient Compressor model --- .../components/compressors/ten_coefficient.py | 24 ++++++++++++------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index cd09349..9adf647 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -1,4 +1,4 @@ -import warnings +import logging from abc import ABC import numpy as np import pandas as pd @@ -6,6 +6,8 @@ from vclibpy.components.compressors.compressor import Compressor from vclibpy.datamodels import Inputs +logger = logging.getLogger(__name__) + def calc_ten_coefficients(T_eva, T_con, coef_list): """ @@ -188,14 +190,14 @@ def get_lambda_h(self, inputs: Inputs): p_outlet = self.get_p_outlet() n_abs = self.get_n_absolute(inputs.n) - T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] + T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T - 273.15 # [°C] if round((self.state_inlet.T - T_eva - 273.15), 2) != round(self.T_sh, 2): - warnings.warn("The superheating of the given state is not " - "equal to the superheating of the datasheet. " - "State1.T_sh= " + str(round((self.state_inlet.T - T_eva - 273.15), 2)) + - ". Datasheet.T_sh = " + str(self.T_sh)) + logger.warning("The superheating of the given state is not " + "equal to the superheating of the datasheet. " + "State1.T_sh= %s. Datasheet.T_sh = %s", + round((self.state_inlet.T - T_eva - 273.15), 2), self.T_sh) # The datasheet has a given superheating temperature which can # vary from the superheating of the real state 1 # which is given by the user. @@ -222,8 +224,11 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): T_con, state_inlet_datasheet, m_flow, capacity, p_el = self._calculate_values( p_2=p_outlet, inputs=inputs ) + if self.T_sc != 0: + h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] + else: + h3 = self.med_prop.calc_state("PQ", p_outlet, 0).h # [J/kg] - h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] h2s = self.med_prop.calc_state("PS", p_outlet, state_inlet_datasheet.s).h # [J/kg] if self._capacity_definition == "heating": @@ -255,8 +260,11 @@ def get_eta_mech(self, inputs: Inputs): T_con, state_inlet_datasheet, m_flow, capacity, p_el = self._calculate_values( p_2=p_outlet, inputs=inputs ) + if self.T_sc != 0: + h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] + else: + h3 = self.med_prop.calc_state("PQ", p_outlet, 0).h # [J/kg] - h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] h2 = h3 + capacity / m_flow # [J/kg] eta_mech = (m_flow * (h2 - state_inlet_datasheet.h)) / p_el From dc8d478d21c9e603e4ba78873ea18ff9610c8036 Mon Sep 17 00:00:00 2001 From: Julius Klupp Date: Wed, 19 Jun 2024 10:21:04 +0200 Subject: [PATCH 19/36] further Bugfixes on TenCoefficient Compressor model --- vclibpy/components/compressors/ten_coefficient.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 9adf647..6820e98 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -198,12 +198,17 @@ def get_lambda_h(self, inputs: Inputs): "equal to the superheating of the datasheet. " "State1.T_sh= %s. Datasheet.T_sh = %s", round((self.state_inlet.T - T_eva - 273.15), 2), self.T_sh) + # The datasheet has a given superheating temperature which can # vary from the superheating of the real state 1 # which is given by the user. # Thus a new self.state_inlet_datasheet has to # be defined for all further calculations - state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) + + if self.T_sh != 0: + state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) + else: + state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 # [kg/s] @@ -284,7 +289,10 @@ def _calculate_values(self, p_2: float, inputs: Inputs): T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] T_con = self.med_prop.calc_state("PQ", p_2, 0).T - 273.15 # [°C] - state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) + if self.T_sh != 0: + state_inlet_datasheet = self.med_prop.calc_state("PT", self.state_inlet.p, T_eva + 273.15 + self.T_sh) + else: + state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 # [kg/s] capacity = self.get_parameter(T_eva, T_con, inputs.n, "capacity") # [W] From f58bf8a3f8ed2ba86de875ba01432c8f16b57191 Mon Sep 17 00:00:00 2001 From: Julius Klupp Date: Wed, 26 Jun 2024 09:17:05 +0200 Subject: [PATCH 20/36] add scaling factor to 10C COmpressor model --- .../components/compressors/ten_coefficient.py | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 6820e98..18881d0 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -167,7 +167,7 @@ class TenCoefficientCompressor(BaseTenCoefficientCompressor): sheet_name (str, optional): Name of the sheet in the datasheet. Defaults to None. """ - def __init__(self, N_max, V_h, T_sc, T_sh, capacity_definition, assumed_eta_mech, datasheet, **kwargs): + def __init__(self, N_max, V_h, T_sc, T_sh, capacity_definition, assumed_eta_mech, datasheet, scaling_factor, **kwargs): super().__init__(N_max=N_max, V_h=V_h, datasheet=datasheet, **kwargs) self.T_sc = T_sc self.T_sh = T_sh @@ -176,6 +176,7 @@ def __init__(self, N_max, V_h, T_sc, T_sh, capacity_definition, assumed_eta_mech self._capacity_definition = capacity_definition self.assumed_eta_mech = assumed_eta_mech self.datasheet = datasheet + self.scaling_factor = scaling_factor def get_lambda_h(self, inputs: Inputs): """ @@ -193,11 +194,11 @@ def get_lambda_h(self, inputs: Inputs): T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T - 273.15 # [°C] - if round((self.state_inlet.T - T_eva - 273.15), 2) != round(self.T_sh, 2): - logger.warning("The superheating of the given state is not " - "equal to the superheating of the datasheet. " - "State1.T_sh= %s. Datasheet.T_sh = %s", - round((self.state_inlet.T - T_eva - 273.15), 2), self.T_sh) + # if round((self.state_inlet.T - T_eva - 273.15), 2) != round(self.T_sh, 2): + # logger.warning("The superheating of the given state is not " + # "equal to the superheating of the datasheet. " + # "State1.T_sh= %s. Datasheet.T_sh = %s", + # round((self.state_inlet.T - T_eva - 273.15), 2), self.T_sh) # The datasheet has a given superheating temperature which can # vary from the superheating of the real state 1 @@ -210,7 +211,7 @@ def get_lambda_h(self, inputs: Inputs): else: state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) - m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 # [kg/s] + m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] # TODO? lambda_h = m_flow / (n_abs * state_inlet_datasheet.d * self.V_h) return lambda_h @@ -294,9 +295,9 @@ def _calculate_values(self, p_2: float, inputs: Inputs): else: state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) - m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 # [kg/s] - capacity = self.get_parameter(T_eva, T_con, inputs.n, "capacity") # [W] - p_el = self.get_parameter(T_eva, T_con, inputs.n, "input_power") # [W] + m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] + capacity = self.get_parameter(T_eva, T_con, inputs.n, "capacity") * self.scaling_factor # [W] + p_el = self.get_parameter(T_eva, T_con, inputs.n, "input_power") * self.scaling_factor # [W] return T_con, state_inlet_datasheet, m_flow, capacity, p_el From 4a91cf1b72e0f3978d57433163fac7a66743c69d Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Wed, 18 Sep 2024 07:40:04 +0200 Subject: [PATCH 21/36] restructure validation repo and function for paper --- tests/test_regression.py | 3 +- .../heat_exchangers/heat_exchanger.py | 1 + vclibpy/utils/nominal_design.py | 41 ++++++++----------- vclibpy/utils/sdf_.py | 4 +- 4 files changed, 22 insertions(+), 27 deletions(-) diff --git a/tests/test_regression.py b/tests/test_regression.py index 3501e0e..9ab0c1f 100644 --- a/tests/test_regression.py +++ b/tests/test_regression.py @@ -155,6 +155,7 @@ def _regression_of_examples(self, flowsheet, fluid): m_flow_eva=0.9, dT_eva_superheating=5, dT_con_subcooling=0, + **kwargs ) path_csv_regression = pathlib.Path(__file__).parent.joinpath( "regression_data", "reference_results", f"{flowsheet}_{fluid}.csv" @@ -171,7 +172,7 @@ def test_vi_ps_propane(self): self._regression_of_examples("VaporInjectionPhaseSeparator", "Propane") def test_evi_propane(self): - self.skipTest("EVI works locally, only CI fails.") + #self.skipTest("EVI works locally, only CI fails.") self._regression_of_examples("VaporInjectionEconomizer", "Propane") @unittest.skip("not implemented") diff --git a/vclibpy/components/heat_exchangers/heat_exchanger.py b/vclibpy/components/heat_exchangers/heat_exchanger.py index 23e8fbe..e56226a 100644 --- a/vclibpy/components/heat_exchangers/heat_exchanger.py +++ b/vclibpy/components/heat_exchangers/heat_exchanger.py @@ -76,6 +76,7 @@ def start_secondary_med_prop(self): def terminate_secondary_med_prop(self): if self.med_prop_sec is not None: self.med_prop_sec.terminate() + self.med_prop_sec = None @abc.abstractmethod def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> Tuple[float, float]: diff --git a/vclibpy/utils/nominal_design.py b/vclibpy/utils/nominal_design.py index 370751b..6a373ab 100644 --- a/vclibpy/utils/nominal_design.py +++ b/vclibpy/utils/nominal_design.py @@ -45,34 +45,27 @@ def nominal_hp_design( m_flow_con_start = kwargs.get("m_flow_con_start", 0.2) m_flow_eva_start = kwargs.get("m_flow_eva_start", 1) accuracy = kwargs.get("accuracy", 0.001) - calculate_m_flow = dT_eva is not None and dT_con is not None - if calculate_m_flow: + + # We have to iterate to match the m_flows to the Q_cons: + m_flow_eva_next = m_flow_eva_start + m_flow_con_next = m_flow_con_start + while True: + # Set values + m_flow_eva = m_flow_eva_next + m_flow_con = m_flow_con_next + inputs.set("m_flow_con", m_flow_con) + inputs.set("m_flow_eva", m_flow_eva) # Get nominal value: fs_state = heat_pump.calc_steady_state(fluid=fluid, inputs=inputs) if fs_state is None: raise ValueError("Given configuration is infeasible at nominal point.") - - else: - # We have to iterate to match the m_flows to the Q_cons: - m_flow_eva_next = m_flow_eva_start - m_flow_con_next = m_flow_con_start - while True: - # Set values - m_flow_eva = m_flow_eva_next - m_flow_con = m_flow_con_next - inputs.set("m_flow_con", m_flow_con) - inputs.set("m_flow_eva", m_flow_eva) - # Get nominal value: - fs_state = heat_pump.calc_steady_state(fluid=fluid, inputs=inputs) - if fs_state is None: - raise ValueError("Given configuration is infeasible at nominal point.") - cp_eva = heat_pump.evaporator._secondary_cp - cp_con = heat_pump.condenser._secondary_cp - m_flow_con_next = fs_state.get("Q_con") / (dT_con * cp_con) - m_flow_eva_next = (fs_state.get("Q_con") * (1 - 1 / fs_state.get("COP"))) / (dT_eva * cp_eva) - # Check convergence: - if abs(m_flow_eva_next - m_flow_eva) < accuracy and abs(m_flow_con-m_flow_con_next) < accuracy: - break + cp_eva = heat_pump.evaporator._secondary_cp + cp_con = heat_pump.condenser._secondary_cp + m_flow_con_next = fs_state.get("Q_con").value / (dT_con * cp_con) + m_flow_eva_next = (fs_state.get("Q_con").value * (1 - 1 / fs_state.get("COP").value)) / (dT_eva * cp_eva) + # Check convergence: + if abs(m_flow_eva_next - m_flow_eva) < accuracy and abs(m_flow_con-m_flow_con_next) < accuracy: + break nominal_design_info = { **inputs.convert_to_str_value_format(with_unit_and_description=False), diff --git a/vclibpy/utils/sdf_.py b/vclibpy/utils/sdf_.py index 1b14d70..470dcfc 100644 --- a/vclibpy/utils/sdf_.py +++ b/vclibpy/utils/sdf_.py @@ -158,6 +158,6 @@ def _unpack_nd_data(data): if __name__ == '__main__': sdf_to_csv( - filepath=pathlib.Path(r"D:\00_temp\calibration_jmo\Optihorst_3D_vclibpy.sdf"), - save_path=pathlib.Path(r"D:\04_git\vclibpy\tests\regression_data") + filepath=pathlib.Path(r"D:\VCLibMap_old.sdf"), + save_path=pathlib.Path(r"D:\00_temp") ) From 2d5720d7528fcd9cd741b7c17c2f77eb059c5136 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Thu, 19 Sep 2024 09:08:00 +0200 Subject: [PATCH 22/36] fix mp error and make calc_multiple_states run in parallel --- vclibpy/flowsheets/base.py | 5 ++++- vclibpy/utils/automation.py | 35 +++++++++++++++++++++++------------ 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index cb046d7..163f29d 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -65,9 +65,12 @@ def setup_new_fluid(self, fluid): self.fluid = fluid def terminate(self): - self.med_prop.terminate() + if self.med_prop is not None: + self.med_prop.terminate() + self.med_prop = None for component in self.get_all_components(): component.terminate_secondary_med_prop() + component.med_prop = None def get_all_components(self) -> List[BaseComponent]: return [self.condenser, self.evaporator] diff --git a/vclibpy/utils/automation.py b/vclibpy/utils/automation.py index 3826ca8..a665421 100644 --- a/vclibpy/utils/automation.py +++ b/vclibpy/utils/automation.py @@ -19,6 +19,8 @@ def calc_multiple_states( save_path: pathlib.Path, heat_pump: BaseCycle, inputs: List[Inputs], + use_multiprocessing: bool = False, + raise_errors: bool = False, **kwargs): """ Function to calculate the flowsheet states for all given inputs. @@ -28,20 +30,27 @@ def calc_multiple_states( save_path (pathlib.Path): Location where to save the results as xlsx. heat_pump (BaseCycle): A valid flowsheet inputs (List[Inputs]): A list with all inputs to simulate + use_multiprocessing (bool): True to use all cores, default no multiprocessing **kwargs: Solver settings for the flowsheet """ rel_infos = [] - for i, single_inputs in enumerate(inputs): - fs_state = None - logger.info(f"Running combination {i+1}/{len(inputs)}.") - try: - fs_state = heat_pump.calc_steady_state(inputs=single_inputs, - **kwargs) - except Exception as e: - # Avoid loss of data if un-excepted errors occur. - logger.error(f"An error occurred: {e}") - if fs_state is None: - fs_state = FlowsheetState() + fs_states = [] + i = 0 + if use_multiprocessing: + mp_inputs = [[heat_pump, inputs_, kwargs, raise_errors] for inputs_ in inputs] + pool = multiprocessing.Pool(processes=min(multiprocessing.cpu_count(), len(inputs))) + for fs_state in pool.imap(_calc_single_hp_state, mp_inputs): + fs_states.append(fs_state) + i += 1 + logger.info(f"Ran {i} of {len(inputs)} points") + else: + for inputs_ in inputs: + fs_state = _calc_single_hp_state([heat_pump, inputs_, kwargs, raise_errors]) + fs_states.append(fs_state) + i += 1 + logger.info(f"Ran {i} of {len(inputs)} points") + + for fs_state, single_inputs in zip(fs_states, inputs): hp_state_dic = { **single_inputs.convert_to_str_value_format(with_unit_and_description=True), **fs_state.convert_to_str_value_format(with_unit_and_description=True) @@ -49,7 +58,9 @@ def calc_multiple_states( rel_infos.append(hp_state_dic) df = pd.DataFrame(rel_infos) df.index.name = "State Number" - df.to_excel(save_path.joinpath(f"{heat_pump}_{heat_pump.fluid}.xlsx"), sheet_name="HP_states", float_format="%.5f") + if os.path.isdir(save_path): + save_path = save_path.joinpath(f"{heat_pump}_{heat_pump.fluid}.xlsx") + df.to_excel(save_path, sheet_name="HP_states", float_format="%.5f") def full_factorial_map_generation( From 991b54f8a9e56ad3507f3d3bcf8e1debfffae398 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Thu, 19 Sep 2024 21:10:19 +0200 Subject: [PATCH 23/36] add kwarg --- vclibpy/utils/automation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/vclibpy/utils/automation.py b/vclibpy/utils/automation.py index a665421..c099857 100644 --- a/vclibpy/utils/automation.py +++ b/vclibpy/utils/automation.py @@ -21,6 +21,7 @@ def calc_multiple_states( inputs: List[Inputs], use_multiprocessing: bool = False, raise_errors: bool = False, + with_unit_and_description: bool = True, **kwargs): """ Function to calculate the flowsheet states for all given inputs. @@ -52,8 +53,8 @@ def calc_multiple_states( for fs_state, single_inputs in zip(fs_states, inputs): hp_state_dic = { - **single_inputs.convert_to_str_value_format(with_unit_and_description=True), - **fs_state.convert_to_str_value_format(with_unit_and_description=True) + **single_inputs.convert_to_str_value_format(with_unit_and_description), + **fs_state.convert_to_str_value_format(with_unit_and_description) } rel_infos.append(hp_state_dic) df = pd.DataFrame(rel_infos) From 1f81baefeb12f890a860d32439d6416d85cb37c0 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Fri, 20 Sep 2024 14:41:34 +0200 Subject: [PATCH 24/36] fixes to get it running --- requirements.txt | 2 +- vclibpy/flowsheets/base.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index 1692d98..5d3fc55 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy>=1.20.0 +numpy<=2.0 matplotlib>=3.3.4 coolprop>=6.4.1 sdf>=0.3.5 diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index 163f29d..c073873 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -180,7 +180,7 @@ def calc_steady_state(self, inputs: Inputs, fluid: str = None, **kwargs): error_eva_history = [] dT_eva_history = [] dT_con_history = [] - error_con, dT_min_eva, dT_min_con, error_eva = np.NAN, np.NAN, np.NAN, np.NAN + error_con, dT_min_eva, dT_min_con, error_eva = np.nan, np.nan, np.nan, np.nan plot_last = -100 if use_quick_solver: From 4487dbb8d148393a73b0cd638a6da1415fdf0f5e Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Thu, 21 Nov 2024 14:44:21 +0100 Subject: [PATCH 25/36] enable option for sdf path saving --- vclibpy/utils/automation.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/vclibpy/utils/automation.py b/vclibpy/utils/automation.py index 2099685..1c0b068 100644 --- a/vclibpy/utils/automation.py +++ b/vclibpy/utils/automation.py @@ -79,6 +79,7 @@ def full_factorial_map_generation( use_multiprocessing: bool = False, save_plots: bool = False, raise_errors: bool = False, + save_sdf: bool = True, **kwargs ) -> (pathlib.Path, pathlib.Path): """ @@ -115,10 +116,12 @@ def full_factorial_map_generation( only relevant if `use_condenser_inlet=False` use_multiprocessing: True to use multiprocessing. May speed up the calculation. Default is False - save_plots: + save_plots (bool): True to save plots of each steady state point. Default is False - raise_errors: + raise_errors (bool): True to raise errors if they occur. + save_sdf (bool): + = False to not save sdf files. Default is True **kwargs: Solver settings for the flowsheet Returns: @@ -198,6 +201,11 @@ def full_factorial_map_generation( save_path_csv, sep=";" ) + # Terminate heat pump med-props: + heat_pump.terminate() + if not save_sdf: + return save_path_csv + for fs_state, idx_triple in zip(fs_states, idx_for_access_later): i_n, i_T_con, i_T_eva_in = idx_triple for variable_name, variable in fs_state.get_variables().items(): @@ -240,10 +248,7 @@ def full_factorial_map_generation( heat_pump.fluid: (_scales, _nd_data, _parameters) } } - # utils.save_to_sdf(data=sdf_data, save_path=save_path_sdf) - - # Terminate heat pump med-props: - heat_pump.terminate() + utils.save_to_sdf(data=sdf_data, save_path=save_path_sdf) return save_path_sdf, save_path_csv From 6ed3a796a0f9cc0805ea4fffec06aaa2f1417804 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Thu, 21 Nov 2024 14:48:11 +0100 Subject: [PATCH 26/36] Add custom eta_mech option, indent --- .../components/compressors/ten_coefficient.py | 38 +++++++++++++------ 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 18881d0..0ecf331 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -1,8 +1,9 @@ import logging from abc import ABC +from typing import Union + import numpy as np import pandas as pd - from vclibpy.components.compressors.compressor import Compressor from vclibpy.datamodels import Inputs @@ -153,7 +154,8 @@ class TenCoefficientCompressor(BaseTenCoefficientCompressor): T_sc (float): Subcooling according to datasheet in K. T_sh (float): Superheating according to datasheet in K. capacity_definition (str): Definition of "capacity" in the datasheet. "cooling" or "heating". - assumed_eta_mech (float): Assumed mechanical efficiency of the compressor (only needed if cooling). + assumed_eta_mech (float, callable): Assumed mechanical efficiency of the compressor (only needed if cooling). + If you pass a funtion, it must have this signature `eta_mech(self, p_outlet, inputs)` datasheet (str): Path of the modified datasheet. **kwargs: parameter_names (dict, optional): @@ -167,14 +169,25 @@ class TenCoefficientCompressor(BaseTenCoefficientCompressor): sheet_name (str, optional): Name of the sheet in the datasheet. Defaults to None. """ - def __init__(self, N_max, V_h, T_sc, T_sh, capacity_definition, assumed_eta_mech, datasheet, scaling_factor, **kwargs): + def __init__( + self, + N_max, V_h, + T_sc, T_sh, capacity_definition, + assumed_eta_mech: Union[float, callable], + datasheet, + scaling_factor, + **kwargs + ): super().__init__(N_max=N_max, V_h=V_h, datasheet=datasheet, **kwargs) self.T_sc = T_sc self.T_sh = T_sh if capacity_definition not in ["cooling", "heating"]: raise ValueError("capacity_definition has to be either 'heating' or 'cooling'") self._capacity_definition = capacity_definition - self.assumed_eta_mech = assumed_eta_mech + if isinstance(assumed_eta_mech, (float, int)): + self.assumed_eta_mech = lambda self, p_outlet, inputs: assumed_eta_mech + else: + self.assumed_eta_mech = assumed_eta_mech self.datasheet = datasheet self.scaling_factor = scaling_factor @@ -191,7 +204,7 @@ def get_lambda_h(self, inputs: Inputs): p_outlet = self.get_p_outlet() n_abs = self.get_n_absolute(inputs.n) - T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] + T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T - 273.15 # [°C] # if round((self.state_inlet.T - T_eva - 273.15), 2) != round(self.T_sh, 2): @@ -211,7 +224,7 @@ def get_lambda_h(self, inputs: Inputs): else: state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) - m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] # TODO? + m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] # TODO? lambda_h = m_flow / (n_abs * state_inlet_datasheet.d * self.V_h) return lambda_h @@ -240,7 +253,10 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): if self._capacity_definition == "heating": h2 = h3 + capacity / m_flow # [J/kg] else: - h2 = h3 + (capacity + p_el * self.assumed_eta_mech) / m_flow # [J/kg] + h2 = h3 + ( + capacity + + p_el * self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) + ) / m_flow # [J/kg] if h2s > h2: raise ValueError("The calculated eta_s is above 1. You probably chose the wrong capacity_definition") @@ -261,7 +277,7 @@ def get_eta_mech(self, inputs: Inputs): p_outlet = self.get_p_outlet() if self._capacity_definition == "cooling": - return self.assumed_eta_mech + return self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) # Else heating T_con, state_inlet_datasheet, m_flow, capacity, p_el = self._calculate_values( p_2=p_outlet, inputs=inputs @@ -295,9 +311,9 @@ def _calculate_values(self, p_2: float, inputs: Inputs): else: state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) - m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] - capacity = self.get_parameter(T_eva, T_con, inputs.n, "capacity") * self.scaling_factor # [W] - p_el = self.get_parameter(T_eva, T_con, inputs.n, "input_power") * self.scaling_factor # [W] + m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] + capacity = self.get_parameter(T_eva, T_con, inputs.n, "capacity") * self.scaling_factor # [W] + p_el = self.get_parameter(T_eva, T_con, inputs.n, "input_power") * self.scaling_factor # [W] return T_con, state_inlet_datasheet, m_flow, capacity, p_el From c35630117785d98f937a8eb3a48dbfea3e3d5eaa Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Thu, 21 Nov 2024 14:48:26 +0100 Subject: [PATCH 27/36] add area-error percentage --- vclibpy/components/heat_exchangers/moving_boundary_Tm.py | 1 + 1 file changed, 1 insertion(+) diff --git a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py index 69c301e..d0db543 100644 --- a/vclibpy/components/heat_exchangers/moving_boundary_Tm.py +++ b/vclibpy/components/heat_exchangers/moving_boundary_Tm.py @@ -129,6 +129,7 @@ def calc(self, inputs: Inputs, fs_state: FlowsheetState) -> (float, float): fs_state.set(name="A_con_sh", value=A_sh, unit="m2", description="Area for superheat heat exchange in condenser") fs_state.set(name="A_con_lat", value=A_lat, unit="m2", description="Area for latent heat exchange in condenser") fs_state.set(name="A_con_sc", value=A_sc, unit="m2", description="Area for subcooling heat exchange in condenser") + fs_state.set(name="error_A", value=error_A, unit="%", description="Area-percentage error for heat exchange in condenser") return error, min(dT_min_in, dT_min_LatSH, From 3b81a2653fa8678eebe8be0a218943e60a1f522f Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Thu, 21 Nov 2024 15:50:56 +0100 Subject: [PATCH 28/36] Only warn and not raise an error if eta_mech assumption is bad --- vclibpy/components/compressors/ten_coefficient.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 0ecf331..877ea85 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -258,10 +258,12 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): p_el * self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) ) / m_flow # [J/kg] - if h2s > h2: - raise ValueError("The calculated eta_s is above 1. You probably chose the wrong capacity_definition") - eta_s = (h2s - state_inlet_datasheet.h) / (h2 - state_inlet_datasheet.h) + if eta_s > 0.8: + logger.warning( + "Calculated eta_is is %s, which is higher than typical maximal values of up to, e.g., 80 %." + "You either chose the wrong capacity_definition, or your assumed eta_mech is also not realistic." + ) return eta_s def get_eta_mech(self, inputs: Inputs): From 30db7d8e273df6d8d352dcbaf4b18762fb1c49c8 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Thu, 21 Nov 2024 15:53:13 +0100 Subject: [PATCH 29/36] add eta_glob to flowsheet --- vclibpy/flowsheets/base.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/vclibpy/flowsheets/base.py b/vclibpy/flowsheets/base.py index c073873..fcdf949 100644 --- a/vclibpy/flowsheets/base.py +++ b/vclibpy/flowsheets/base.py @@ -501,6 +501,10 @@ def calculate_outputs_for_valid_pressures( name="dT_min_con", value=dT_min_con, unit="K", description="Condenser pinch temperature" ) + fs_state.set( + name="eta_glob", value=fs_state.get("eta_is").value * fs_state.get("eta_mech").value, + unit="%", description="Global compressor efficiency" + ) if save_path_plots is not None: input_name = inputs.get_name() self.plot_cycle(save_path=save_path_plots.joinpath(f"{input_name}_final_result.png"), inputs=inputs) From cc885eae854898308c569c0183199366fa390fa5 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Fri, 22 Nov 2024 16:30:54 +0100 Subject: [PATCH 30/36] use eta_is and don't save as xlsx to avoid need for xls writer. Use compressor instance directly --- .../components/compressors/ten_coefficient.py | 20 +++-- .../ten_coefficient_compressor_reqression.py | 87 +++++++------------ 2 files changed, 40 insertions(+), 67 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 877ea85..bcb786f 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -58,7 +58,7 @@ class BaseTenCoefficientCompressor(Compressor, ABC): "m_flow": "Flow Rate(kg/h)", "capacity": "Capacity(W)", "input_power": "Input Power(W)", - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumentric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } @@ -86,7 +86,7 @@ def __init__(self, N_max, V_h, datasheet, **kwargs): "m_flow": "Flow Rate(kg/h)", "capacity": "Capacity(W)", "input_power": "Input Power(W)", - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumentric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } @@ -258,13 +258,15 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): p_el * self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) ) / m_flow # [J/kg] - eta_s = (h2s - state_inlet_datasheet.h) / (h2 - state_inlet_datasheet.h) - if eta_s > 0.8: + eta_is = (h2s - state_inlet_datasheet.h) / (h2 - state_inlet_datasheet.h) + if eta_is > 0.8: logger.warning( - "Calculated eta_is is %s, which is higher than typical maximal values of up to, e.g., 80 %." - "You either chose the wrong capacity_definition, or your assumed eta_mech is also not realistic." + f"Calculated eta_is is {eta_is * 100} %, which is higher than " + f"typical maximal values of up to, e.g., 80 %." + "You either chose the wrong capacity_definition, " + "or your assumed eta_mech is also not realistic.", ) - return eta_s + return eta_is def get_eta_mech(self, inputs: Inputs): """ @@ -346,7 +348,7 @@ class DataSheetCompressor(BaseTenCoefficientCompressor): Dictionary to match internal parameter names (keys) to the names used in the table values. Default { - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumetric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } @@ -384,7 +386,7 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): """ T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 0).T T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T - return self.get_parameter(T_eva, T_con, inputs.n, "eta_s") + return self.get_parameter(T_eva, T_con, inputs.n, "eta_is") def get_eta_mech(self, inputs: Inputs): """ diff --git a/vclibpy/utils/ten_coefficient_compressor_reqression.py b/vclibpy/utils/ten_coefficient_compressor_reqression.py index a3b484e..5a3f908 100644 --- a/vclibpy/utils/ten_coefficient_compressor_reqression.py +++ b/vclibpy/utils/ten_coefficient_compressor_reqression.py @@ -1,6 +1,7 @@ -from typing import List +from typing import List, Union import csv import pandas as pd +from pathlib import Path from vclibpy.components.compressors import TenCoefficientCompressor from vclibpy import media, Inputs @@ -8,72 +9,58 @@ try: from sklearn.linear_model import LinearRegression from sklearn.preprocessing import PolynomialFeatures - from xlsxwriter.workbook import Workbook except ImportError as err: - raise ImportError("You have to install xlsxwriter and " - "sklearn to use this regression tool") + raise ImportError("You have to install sklearn to use this regression tool") def create_regression_data( - variables: List[str], - T_con: List[float], T_eva: List[float], n: List[float], - T_sh: float, T_sc: float, n_max: int, V_h: float, fluid: str, datasheet: str, - capacity_definition: str, assumed_eta_mech: int, - folder_path: str): + T_con: List[float], + T_eva: List[float], + n: List[float], + compressor: TenCoefficientCompressor, + save_path: Union[Path, str], + fluid: str = None, +): """ Performs multidimensional linear regression to create compressor learning data. Args: - variables: (List[str]): - Variable names to create regressions for. - Options are: eta_s, lambda_h, and eta_mech T_con (List[float]): Condensing temperatures in K T_eva (List[float]): Evaporation temperatures in K n (List[float]): Compressor speeds from 0 to 1 - T_sh (float): Superheat temperature. - T_sc (float): Subcooling temperature. - n_max (int): Maximum compressor speed. - V_h (float): Compressor volume. - fluid (str): Type of refrigerant. - datasheet (str): Path to the modified datasheet. - capacity_definition (str): Definition of compressor capacity (e.g., "cooling"). - assumed_eta_mech (int): Assumed mechanical efficiency. - folder_path (str): Path to the folder where the newly created table will be saved. + save_path (str): Path to the folder where the newly created table will be saved. + compressor (TenCoefficientCompressor): + Instance to create regression for. + If med_prop is not set, the given fluid will be used. + fluid (str): Type of refrigerant. Only used if not already specified in compressor. Returns: List[float]: A list containing the regression parameters [P0, P1, ..., P9]. Raises: ValueError: If any specified variable column is not present in the DataFrame. - - Example: - >>> create_regression_data(11.1, 8.3, 120, 20.9e-6, "PROPANE", "path/to/datasheet.xlsx", - ... "cooling", 1, 6, 5, 5, 5, 303.15, 243.15, "path/to/save", 3) - [intercept, P1, P2, P3, P4, P5, P6, P7, P8, P9] """ - # create RefProp, fluid & compressor instance - med_prop = media.CoolProp(fluid) - - compressor = TenCoefficientCompressor( - N_max=n_max, V_h=V_h, T_sc=T_sc, T_sh=T_sh, - capacity_definition=capacity_definition, - assumed_eta_mech=assumed_eta_mech, - datasheet=datasheet - ) + if compressor.med_prop is None: + med_prop_class, med_prop_kwargs = media.get_global_med_prop_and_kwargs() + med_prop = med_prop_class(fluid_name=fluid, **med_prop_kwargs) + compressor.med_prop = med_prop + else: + med_prop = compressor.med_prop - compressor.med_prop = med_prop keywords = { - "eta_s": "Isentropic Efficiency(-)", + "eta_is": "Isentropic Efficiency(-)", "lambda_h": "Volumentric Efficiency(-)", "eta_mech": "Mechanical Efficiency(-)" } + variables = ["eta_is", "eta_mech", "lambda_h"] + tuples_for_cols = [("", "n")] for _variable in variables: for _n in n: tuples_for_cols.append((keywords[_variable], compressor.get_n_absolute(_n))) # tuples_for_cols: - # eta_s eta_s eta_s lambda_h lambda_h lambda_h eta_mech ... + # eta_is eta_is eta_is lambda_h lambda_h lambda_h eta_mech ... # n 30 60 90 30 60 90 30 ... cols = pd.MultiIndex.from_tuples(tuples_for_cols) final_df = pd.DataFrame( @@ -83,7 +70,7 @@ def create_regression_data( # final_df: column names are tuples (tuples_for_cols). # First column is filled with P1, P2, ... - # for-loop for multiple types(eta_s, eta_mech, etc) + # for-loop for multiple types(eta_is, eta_mech, etc) for m, _variable in enumerate(variables): for k, _n in enumerate(n): # for-loop for multiple rotation speeds T_eva_list = [] @@ -95,7 +82,7 @@ def create_regression_data( for j in range(len(T_con)): if T_eva[i] < T_con[j]: p1 = med_prop.calc_state("TQ", T_eva[i], 1).p - state_1 = med_prop.calc_state("PT", p1, (T_eva[i] + T_sh)) + state_1 = med_prop.calc_state("PT", p1, (T_eva[i] + compressor.T_sh)) compressor.state_inlet = state_1 p2 = med_prop.calc_state("TQ", T_con[j], 1).p # The enthalpy and entropy of the outlet @@ -107,7 +94,7 @@ def create_regression_data( T_con_list.append(T_con[j]) inputs = Inputs(n=_n) - if _variable == "eta_s": + if _variable == "eta_is": result_list.append(compressor.get_eta_isentropic( p_outlet=p2, inputs=inputs) ) @@ -123,23 +110,7 @@ def create_regression_data( ) final_df[cols[m * len(n) + k + 1]] = create_regression_parameters(df, _variable) - - # dataframes with a double column header can't be saved as a - # .xlsx yet, if index = False (NotImplementedError). - # .xlsx format is necessary, because TenCoefficientCompressor.get_parameter() - # expects a .xlsx as an input - # --> workaround: save the dataframe as a .csv, read it again and save it as a .xlsx - # TODO: Revert once this feature is in pandas. - final_df.to_csv(folder_path + r"\regressions.csv", index=False) - - workbook = Workbook(folder_path + r"\regressions.xlsx") - worksheet = workbook.add_worksheet() - with open(folder_path + r"\regressions.csv", 'rt', encoding='utf8') as f: - reader = csv.reader(f) - for r, row in enumerate(reader): - for c, col in enumerate(row): - worksheet.write(r, c, col) - workbook.close() + final_df.to_csv(save_path, index=False) def create_regression_parameters(df: pd.DataFrame, variable: str): From fb0b5a57a2383571176457f43651020af8366e84 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Sat, 23 Nov 2024 12:02:18 +0100 Subject: [PATCH 31/36] enable csv input and add extrapolation --- .../components/compressors/ten_coefficient.py | 50 ++++++++++++++++--- .../ten_coefficient_compressor_reqression.py | 4 ++ 2 files changed, 46 insertions(+), 8 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index bcb786f..3c8abe7 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -63,23 +63,22 @@ class BaseTenCoefficientCompressor(Compressor, ABC): "eta_mech": "Mechanical Efficiency(-)" } sheet_name (str, optional): Name of the sheet in the datasheet. Defaults to None. + extrapolate (str, optional): + Method to handle extrapolation of data. + Default "hold" means no extrapolation """ def __init__(self, N_max, V_h, datasheet, **kwargs): """ Initialize the BaseTenCoefficientCompressor. - - Args: - N_max (float): Maximal rotations per second of the compressor. - V_h (float): Volume of the compressor in m^3. - datasheet (str): Path of the datasheet file. - parameter_names (dict, optional): Dictionary of parameter names. Defaults to None. - sheet_name (str, optional): Name of the sheet in the datasheet. Defaults to None. """ super().__init__(N_max, V_h) sheet_name = kwargs.get('sheet_name', None) - self.md = pd.read_excel(datasheet, sheet_name=sheet_name) + if str(datasheet).endswith(".xslx"): + self.md = pd.read_excel(datasheet, sheet_name=sheet_name) + else: + self.md = pd.read_csv(datasheet) parameter_names = kwargs.get('parameter_names', None) if parameter_names is None: self.parameter_names = { @@ -92,6 +91,7 @@ def __init__(self, N_max, V_h, datasheet, **kwargs): } else: self.parameter_names = parameter_names + self.extrapolate = kwargs.get("extrapolate", "hold") def get_parameter(self, T_eva, T_con, n, type_): """ @@ -122,6 +122,14 @@ def get_parameter(self, T_eva, T_con, n, type_): return np.interp(self.get_n_absolute(n), n_list, param_list) # linear interpolation + def _interpolate(self): + if self.extrapolate == "hold": + # linear interpolation, no extrapolation + return np.interp(self.get_n_absolute(n), n_list, param_list) + if self.extrapolate == "linear": + return linear_interpolate_extrapolate(self.get_n_absolute(n), n_list, param_list) + raise KeyError(f"Given extrapolate option '{self.extrapolate}' is not supported!") + class TenCoefficientCompressor(BaseTenCoefficientCompressor): """ @@ -402,3 +410,29 @@ def get_eta_mech(self, inputs: Inputs): T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 0).T T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T return self.get_parameter(T_eva, T_con, inputs.n, "eta_mech") + + +def linear_interpolate_extrapolate(x_new, x, y): + """ + Util function for linear 1D extrapolation or interpolation. + Used to avoid scipy requirement. + + x_new: points where to interpolate/extrapolate + x: known x values + y: known y values + """ + y_new = np.interp(x_new, x, y) + + # Handle left extrapolation + left_mask = x_new < x[0] + if np.any(left_mask): + slope = (y[1] - y[0]) / (x[1] - x[0]) + y_new[left_mask] = y[0] + slope * (x_new[left_mask] - x[0]) + + # Handle right extrapolation + right_mask = x_new > x[-1] + if np.any(right_mask): + slope = (y[-1] - y[-2]) / (x[-1] - x[-2]) + y_new[right_mask] = y[-1] + slope * (x_new[right_mask] - x[-1]) + + return y_new diff --git a/vclibpy/utils/ten_coefficient_compressor_reqression.py b/vclibpy/utils/ten_coefficient_compressor_reqression.py index 5a3f908..3722e68 100644 --- a/vclibpy/utils/ten_coefficient_compressor_reqression.py +++ b/vclibpy/utils/ten_coefficient_compressor_reqression.py @@ -110,6 +110,10 @@ def create_regression_data( ) final_df[cols[m * len(n) + k + 1]] = create_regression_parameters(df, _variable) + df_new = final_df.copy() + df_new.columns = df_new.columns.get_level_values(0) # Use only first level + df_new.loc[-1] = final_df.columns.get_level_values(1) + final_df = df_new.sort_index().reset_index(drop=True) final_df.to_csv(save_path, index=False) From 912c3556c277e98d4dab00af6199cff7cc649449 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Sat, 23 Nov 2024 12:11:02 +0100 Subject: [PATCH 32/36] bugfix array --- .../components/compressors/ten_coefficient.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 3c8abe7..5f6b7d6 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -120,14 +120,14 @@ def get_parameter(self, T_eva, T_con, n, type_): n_list.append(coefficients.pop(0)) param_list.append(calc_ten_coefficients(T_eva, T_con, coefficients)) - return np.interp(self.get_n_absolute(n), n_list, param_list) # linear interpolation + return self._interpolate(self.get_n_absolute(n), n_list, param_list) # linear interpolation - def _interpolate(self): + def _interpolate(self, x_new, x, y): if self.extrapolate == "hold": # linear interpolation, no extrapolation - return np.interp(self.get_n_absolute(n), n_list, param_list) + return np.interp(x_new, x, y) if self.extrapolate == "linear": - return linear_interpolate_extrapolate(self.get_n_absolute(n), n_list, param_list) + return linear_interpolate_extrapolate(x_new, x, y) raise KeyError(f"Given extrapolate option '{self.extrapolate}' is not supported!") @@ -425,14 +425,13 @@ def linear_interpolate_extrapolate(x_new, x, y): # Handle left extrapolation left_mask = x_new < x[0] - if np.any(left_mask): + if x_new < x[0]: slope = (y[1] - y[0]) / (x[1] - x[0]) - y_new[left_mask] = y[0] + slope * (x_new[left_mask] - x[0]) + y_new = y[0] + slope * (x_new - x[0]) # Handle right extrapolation - right_mask = x_new > x[-1] - if np.any(right_mask): + if x_new > x[-1]: slope = (y[-1] - y[-2]) / (x[-1] - x[-2]) - y_new[right_mask] = y[-1] + slope * (x_new[right_mask] - x[-1]) + y_new = y[-1] + slope * (x_new - x[-1]) return y_new From c2080ca6d4021f8e0425592f23e06c8c91655d10 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Sat, 23 Nov 2024 12:56:49 +0100 Subject: [PATCH 33/36] fix xlsx --- .../components/compressors/ten_coefficient.py | 22 ++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 5f6b7d6..709d086 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -75,7 +75,7 @@ def __init__(self, N_max, V_h, datasheet, **kwargs): super().__init__(N_max, V_h) sheet_name = kwargs.get('sheet_name', None) - if str(datasheet).endswith(".xslx"): + if str(datasheet).endswith(".xlsx"): self.md = pd.read_excel(datasheet, sheet_name=sheet_name) else: self.md = pd.read_csv(datasheet) @@ -192,10 +192,7 @@ def __init__( if capacity_definition not in ["cooling", "heating"]: raise ValueError("capacity_definition has to be either 'heating' or 'cooling'") self._capacity_definition = capacity_definition - if isinstance(assumed_eta_mech, (float, int)): - self.assumed_eta_mech = lambda self, p_outlet, inputs: assumed_eta_mech - else: - self.assumed_eta_mech = assumed_eta_mech + self.assumed_eta_mech = assumed_eta_mech self.datasheet = datasheet self.scaling_factor = scaling_factor @@ -258,13 +255,15 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): h2s = self.med_prop.calc_state("PS", p_outlet, state_inlet_datasheet.s).h # [J/kg] + if callable(self.assumed_eta_mech): + eta_mech = self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) + else: + eta_mech = self.assumed_eta_mech + if self._capacity_definition == "heating": h2 = h3 + capacity / m_flow # [J/kg] else: - h2 = h3 + ( - capacity + - p_el * self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) - ) / m_flow # [J/kg] + h2 = h3 + (capacity + p_el * eta_mech) / m_flow # [J/kg] eta_is = (h2s - state_inlet_datasheet.h) / (h2 - state_inlet_datasheet.h) if eta_is > 0.8: @@ -289,7 +288,10 @@ def get_eta_mech(self, inputs: Inputs): p_outlet = self.get_p_outlet() if self._capacity_definition == "cooling": - return self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) + if isinstance(self.assumed_eta_mech, callable): + return self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) + return self.assumed_eta_mech + # Else heating T_con, state_inlet_datasheet, m_flow, capacity, p_el = self._calculate_values( p_2=p_outlet, inputs=inputs From 4126fcf2f7842719a8c201ab52b9f4525cc28308 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Sat, 23 Nov 2024 12:59:49 +0100 Subject: [PATCH 34/36] add comment on why lambda is not used --- vclibpy/components/compressors/ten_coefficient.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 709d086..e121a5b 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -192,6 +192,8 @@ def __init__( if capacity_definition not in ["cooling", "heating"]: raise ValueError("capacity_definition has to be either 'heating' or 'cooling'") self._capacity_definition = capacity_definition + # Don't use lambda function to cast float as a function, + # as local functions are not pickable for multiprocessing self.assumed_eta_mech = assumed_eta_mech self.datasheet = datasheet self.scaling_factor = scaling_factor @@ -288,7 +290,7 @@ def get_eta_mech(self, inputs: Inputs): p_outlet = self.get_p_outlet() if self._capacity_definition == "cooling": - if isinstance(self.assumed_eta_mech, callable): + if callable(self.assumed_eta_mech): return self.assumed_eta_mech(self=self, p_outlet=p_outlet, inputs=inputs) return self.assumed_eta_mech From bc3efa9e554693c935d8f1adabc43e89c73053be Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Sun, 24 Nov 2024 18:01:00 +0100 Subject: [PATCH 35/36] dont require T_sc for cooling definition --- .../components/compressors/ten_coefficient.py | 30 +++++++++++-------- 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index e121a5b..6e7a4dc 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -179,14 +179,21 @@ class TenCoefficientCompressor(BaseTenCoefficientCompressor): def __init__( self, - N_max, V_h, - T_sc, T_sh, capacity_definition, - assumed_eta_mech: Union[float, callable], + N_max, + V_h, + T_sh, + capacity_definition, datasheet, - scaling_factor, + T_sc: float = None, + assumed_eta_mech: Union[float, callable] = None, + scaling_factor: float = 1, **kwargs ): super().__init__(N_max=N_max, V_h=V_h, datasheet=datasheet, **kwargs) + if capacity_definition == "cooling" and assumed_eta_mech is None: + raise ValueError("capacity_definition cooling requires an assumption for eta_mech") + if capacity_definition == "heating" and T_sc is None: + raise ValueError("capacity_definition heating requires an assumption for T_sc") self.T_sc = T_sc self.T_sh = T_sh if capacity_definition not in ["cooling", "heating"]: @@ -195,7 +202,6 @@ def __init__( # Don't use lambda function to cast float as a function, # as local functions are not pickable for multiprocessing self.assumed_eta_mech = assumed_eta_mech - self.datasheet = datasheet self.scaling_factor = scaling_factor def get_lambda_h(self, inputs: Inputs): @@ -231,7 +237,7 @@ def get_lambda_h(self, inputs: Inputs): else: state_inlet_datasheet = self.med_prop.calc_state("PQ", self.state_inlet.p, 1) - m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] # TODO? + m_flow = self.get_parameter(T_eva, T_con, inputs.n, "m_flow") / 3600 * self.scaling_factor # [kg/s] lambda_h = m_flow / (n_abs * state_inlet_datasheet.d * self.V_h) return lambda_h @@ -250,10 +256,6 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): T_con, state_inlet_datasheet, m_flow, capacity, p_el = self._calculate_values( p_2=p_outlet, inputs=inputs ) - if self.T_sc != 0: - h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] - else: - h3 = self.med_prop.calc_state("PQ", p_outlet, 0).h # [J/kg] h2s = self.med_prop.calc_state("PS", p_outlet, state_inlet_datasheet.s).h # [J/kg] @@ -263,15 +265,19 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): eta_mech = self.assumed_eta_mech if self._capacity_definition == "heating": + if self.T_sc != 0: + h3 = self.med_prop.calc_state("PT", p_outlet, T_con + 273.15 - self.T_sc).h # [J/kg] + else: + h3 = self.med_prop.calc_state("PQ", p_outlet, 0).h # [J/kg] h2 = h3 + capacity / m_flow # [J/kg] else: - h2 = h3 + (capacity + p_el * eta_mech) / m_flow # [J/kg] + h2 = state_inlet_datasheet.h + (p_el * eta_mech) / m_flow # [J/kg] eta_is = (h2s - state_inlet_datasheet.h) / (h2 - state_inlet_datasheet.h) if eta_is > 0.8: logger.warning( f"Calculated eta_is is {eta_is * 100} %, which is higher than " - f"typical maximal values of up to, e.g., 80 %." + f"typical maximal values of up to, e.g., 80 %. " "You either chose the wrong capacity_definition, " "or your assumed eta_mech is also not realistic.", ) From 41bc1d3550360706addc4ca94bca562b278bce00 Mon Sep 17 00:00:00 2001 From: "fabian.wuellhorst" Date: Mon, 25 Nov 2024 13:59:04 +0100 Subject: [PATCH 36/36] add changelog and increase version --- CHANGELOG.md | 11 +++++++++++ vclibpy/__init__.py | 2 +- vclibpy/components/compressors/ten_coefficient.py | 12 ++++++------ 3 files changed, 18 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 07216be..64c5357 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,2 +1,13 @@ - v0.1.0: - First implementation based on prior work +- v0.2.0: + - Add LMTD MovingBoundary heat exchangers + - Add VDI atlas methods for plate heat exchanger + - Add scaling factor to TenCoefficientCompressor + - Add extrapolation option to TenCoefficientCompressor + - Improve warnings and parameters in TenCoefficientCompressor (eta_mech, T_sc) + - Add T_con_out as possible Input + - Automation is now with multiprocessing + - SDF saving is optional in automation + - Improve nominal design util + - Improve C10 regression generation, add .csv support for TenCoefficientCompressor diff --git a/vclibpy/__init__.py b/vclibpy/__init__.py index 6bb4742..1be17e7 100644 --- a/vclibpy/__init__.py +++ b/vclibpy/__init__.py @@ -4,4 +4,4 @@ from vclibpy.datamodels import FlowsheetState, Inputs -__version__ = "0.1.0" +__version__ = "0.2.0" diff --git a/vclibpy/components/compressors/ten_coefficient.py b/vclibpy/components/compressors/ten_coefficient.py index 6e7a4dc..f6a87c5 100644 --- a/vclibpy/components/compressors/ten_coefficient.py +++ b/vclibpy/components/compressors/ten_coefficient.py @@ -220,11 +220,11 @@ def get_lambda_h(self, inputs: Inputs): T_eva = self.med_prop.calc_state("PQ", self.state_inlet.p, 1).T - 273.15 # [°C] T_con = self.med_prop.calc_state("PQ", p_outlet, 0).T - 273.15 # [°C] - # if round((self.state_inlet.T - T_eva - 273.15), 2) != round(self.T_sh, 2): - # logger.warning("The superheating of the given state is not " - # "equal to the superheating of the datasheet. " - # "State1.T_sh= %s. Datasheet.T_sh = %s", - # round((self.state_inlet.T - T_eva - 273.15), 2), self.T_sh) + if round((self.state_inlet.T - T_eva - 273.15), 2) != round(self.T_sh, 2): + logger.debug("The superheating of the given state is not " + "equal to the superheating of the datasheet. " + "State: T_sh=%s; Datasheet: T_sh=%s", + round((self.state_inlet.T - T_eva - 273.15), 2), self.T_sh) # The datasheet has a given superheating temperature which can # vary from the superheating of the real state 1 @@ -279,7 +279,7 @@ def get_eta_isentropic(self, p_outlet: float, inputs: Inputs): f"Calculated eta_is is {eta_is * 100} %, which is higher than " f"typical maximal values of up to, e.g., 80 %. " "You either chose the wrong capacity_definition, " - "or your assumed eta_mech is also not realistic.", + "or your assumed eta_mech is not realistic.", ) return eta_is