diff --git a/src/tespy/components/heat_exchangers/base.py b/src/tespy/components/heat_exchangers/base.py index a787374ce..d565b912b 100644 --- a/src/tespy/components/heat_exchangers/base.py +++ b/src/tespy/components/heat_exchangers/base.py @@ -465,9 +465,9 @@ def calculate_td_log(self): T_o2 = o2.calc_T() if T_i1 <= T_o2: - T_i1 = T_o2 + 0.01 + T_o2 = T_i1 - 0.1 if T_o1 <= T_i2: - T_o1 = T_i2 + 0.01 + T_o1 = T_i2 + 0.1 ttd_u = T_i1 - T_o2 ttd_l = T_o1 - T_i2 @@ -563,8 +563,17 @@ def kA_char_func(self): """ p1 = self.kA_char1.param p2 = self.kA_char2.param - f1 = self.get_char_expr(p1, **self.kA_char1.char_params) - f2 = self.get_char_expr(p2, **self.kA_char2.char_params) + if self.local_offdesign: + design_value = self._connection_offdesign[self.inl[0].label][p1] + actual_value = getattr(self.inl[0], p1).val_SI + f1 = actual_value / design_value + + design_value = self._connection_offdesign[self.inl[1].label][p2] + actual_value = getattr(self.inl[1], p2).val_SI + f2 = actual_value / design_value + else: + f1 = self.get_char_expr(p1, **self.kA_char1.char_params) + f2 = self.get_char_expr(p2, **self.kA_char2.char_params) fkA1 = self.kA_char1.char_func.evaluate(f1) fkA2 = self.kA_char2.char_func.evaluate(f2) diff --git a/src/tespy/components/heat_exchangers/sectioned.py b/src/tespy/components/heat_exchangers/sectioned.py index 082961a9f..58bdfe7b2 100644 --- a/src/tespy/components/heat_exchangers/sectioned.py +++ b/src/tespy/components/heat_exchangers/sectioned.py @@ -759,9 +759,9 @@ def UA_cecchinato_func(self): secondary_index = 0 m_r = self.inl[refrigerant_index].m - m_ratio_r = m_r.val_SI / m_r.design + m_ratio_r = max(m_r.val_SI / m_r.design, 1e-6) m_sf = self.inl[secondary_index].m - m_ratio_sf = m_sf.val_SI / m_sf.design + m_ratio_sf = max(m_sf.val_SI / m_sf.design, 1e-6) fUA = ( (1 + alpha_ratio * area_ratio) diff --git a/src/tespy/components/nodes/humidity_control.py b/src/tespy/components/nodes/humidity_control.py new file mode 100644 index 000000000..1d63d4279 --- /dev/null +++ b/src/tespy/components/nodes/humidity_control.py @@ -0,0 +1,324 @@ +# -*- coding: utf-8 + +"""Module of class HumidityControl. + + +This file is part of project TESPy (github.com/oemof/tespy). It's copyrighted +by the contributors recorded in the version control history of the file, +available from its original location tespy/components/nodes/humidity_control.py + +SPDX-License-Identifier: MIT +""" + +from tespy.components.component import component_registry +from tespy.components.nodes.base import NodeBase +from tespy.tools.data_containers import ComponentMandatoryConstraints as dc_cmc +from tespy.tools.data_containers import SimpleDataContainer as dc_simple +from tespy.tools.fluid_properties import dT_mix_dph +from tespy.tools.fluid_properties import dT_mix_pdh + +# from tespy.tools.fluid_properties import dT_mix_ph_dfluid + + +@component_registry +class HumidityControl(NodeBase): + r""" + A humidity control handles the water content from a humid air flow. + + **Mandatory Equations** + + - :py:meth:`tespy.components.nodes.base.NodeBase.mass_flow_func` + - :py:meth:`tespy.components.nodes.base.NodeBase.pressure_structure_matrix` + - :py:meth:`tespy.components.nodes.humidity_control.HumidityControl.fluid_func` + - :py:meth:`tespy.components.nodes.humidity_control.HumidityControl.energy_balance_func` + + Inlets/Outlets + + - in1: inlet with humid air + - out1: outlet with humid air + - out2: outlet of liquid water or ice + + Image + + .. image:: /api/_images/Splitter.svg + :alt: flowsheet of the splitter + :align: center + :class: only-light + + .. image:: /api/_images/Splitter_darkmode.svg + :alt: flowsheet of the splitter + :align: center + :class: only-dark + + Note + ---- + Fluid separation requires power and cooling, equations have not been + implemented, yet! + + Parameters + ---------- + label : str + The label of the component. + + design : list + List containing design parameters (stated as String). + + offdesign : list + List containing offdesign parameters (stated as String). + + design_path : str + Path to the components design case. + + local_offdesign : boolean + Treat this component in offdesign mode in a design calculation. + + local_design : boolean + Treat this component in design mode in an offdesign calculation. + + char_warnings : boolean + Ignore warnings on default characteristics usage for this component. + + printout : boolean + Include this component in the network's results printout. + + num_out : float, dict + Number of outlets for this component, default value: 2. + + Example + ------- + The separator is used to split up a single mass flow into a specified + number of different parts at identical pressure and temperature but + different fluid composition. Fluids can be separated from each other. + + >>> from tespy.components import Sink, Source, Separator + >>> from tespy.connections import Connection + >>> from tespy.networks import Network + >>> nw = Network(iterinfo=False) + >>> nw.units.set_defaults(**{ + ... "pressure": "bar", "temperature": "degC" + ... }) + >>> so = Source('source') + >>> si1 = Sink('sink1') + >>> si2 = Sink('sink2') + >>> s = Separator('separator', num_out=2) + >>> inc = Connection(so, 'out1', s, 'in1') + >>> outg1 = Connection(s, 'out1', si1, 'in1') + >>> outg2 = Connection(s, 'out2', si2, 'in1') + >>> nw.add_conns(inc, outg1, outg2) + + An Air (simplified) mass flow of 5 kg/s is split up into two mass flows. + One mass flow of 1 kg/s containing 10 % oxygen and 90 % nitrogen leaves the + separator. It is possible to calculate the fluid composition of the second + mass flow. Specify starting values for the second mass flow fluid + composition for calculation stability. + + >>> inc.set_attr(fluid={'O2': 0.23, 'N2': 0.77}, p=1, T=20, m=5) + >>> outg1.set_attr(fluid={'O2': 0.1, 'N2': 0.9}, m=1) + >>> outg2.set_attr(fluid0={'O2': 0.5, 'N2': 0.5}) + >>> nw.solve('design') + >>> outg2.fluid.val['O2'] + 0.2625 + + In the same way, it is possible to specify one of the fluid components in + the second mass flow instead of the first mass flow. The solver will find + the mass flows matching the desired composition. 65 % of the mass flow + will leave the separator at the second outlet the case of 30 % oxygen + mass fraction for this outlet. + + >>> outg1.set_attr(m=None) + >>> outg2.set_attr(fluid={'O2': 0.3}) + >>> nw.solve('design') + >>> outg2.fluid.val['O2'] + 0.3 + >>> round(outg2.m.val_SI / inc.m.val_SI, 2) + 0.65 + """ + + @staticmethod + def get_parameters(): + return {'num_out': dc_simple(description="number of outlets")} + + def _update_num_eq(self): + self.variable_fluids = set( + [fluid for c in self.inl + self.outl for fluid in c.fluid.is_var] + ) + num_fluid_eq = len(self.variable_fluids) + if num_fluid_eq == 0: + num_fluid_eq = 1 + self.variable_fluids = [list(self.inl[0].fluid.is_set)[0]] + + self.constraints["fluid_constraints"].num_eq = num_fluid_eq + + def get_mandatory_constraints(self): + return { + 'mass_flow_constraints': dc_cmc(**{ + 'num_eq_sets': 1, + 'func': self.mass_flow_func, + 'dependents': self.mass_flow_dependents, + 'description': 'mass balance constraint' + }), + 'fluid_constraints': dc_cmc(**{ + 'num_eq_sets': self.num_o, + 'func': self.fluid_func, + 'deriv': self.fluid_deriv, + 'dependents': self.fluid_dependents, + 'description': 'fluid mass fraction balance constraints' + }), + 'energy_balance_constraints': dc_cmc(**{ + 'num_eq_sets': self.num_o, + 'func': self.energy_balance_func, + 'deriv': self.energy_balance_deriv, + 'dependents': self.energy_balance_dependents, + 'description': 'equal temperature at all outlets constraints' + }), + 'pressure_constraints': dc_cmc(**{ + 'num_eq_sets': self.num_o, + 'structure_matrix': self.pressure_structure_matrix, + 'description': 'pressure equality constraints' + }) + } + + @staticmethod + def inlets(): + return ['in1'] + + def outlets(self): + if self.num_out.is_set: + return [f'out{i + 1}' for i in range(self.num_out.val)] + else: + self.set_attr(num_out=2) + return self.outlets() + + def propagate_wrapper_to_target(self, branch): + branch["components"] += [self] + for outconn in self.outl: + branch["connections"] += [outconn] + outconn.target.propagate_wrapper_to_target(branch) + + def fluid_func(self): + r""" + Calculate the vector of residual values for fluid balance equations. + + Returns + ------- + residual : list + Vector of residual values for component's fluid balance. + + .. math:: + + 0 = \dot{m}_{in} \cdot x_{fl,in} - \dot {m}_{out,j} + \cdot x_{fl,out,j}\\ + \forall fl \in \text{network fluids,} + \; \forall j \in \text{outlets} + """ + i = self.inl[0] + residual = [] + for fluid in self.variable_fluids: + res = i.fluid.val[fluid] * i.m.val_SI + for o in self.outl: + res -= o.fluid.val[fluid] * o.m.val_SI + residual += [res] + return residual + + def fluid_deriv(self, increment_filter, k, dependents=None): + r""" + Calculate partial derivatives of fluid balance. + + Parameters + ---------- + increment_filter : ndarray + Matrix for filtering non-changing variables. + + k : int + Position of derivatives in Jacobian matrix (k-th equation). + """ + i = self.inl[0] + for fluid in self.variable_fluids: + for o in self.outl: + self._partial_derivative(o.m, k, -o.fluid.val[fluid], increment_filter) + if fluid in o.fluid.is_var: + self.jacobian[k, o.fluid.J_col[fluid]] = -o.m.val_SI + + self._partial_derivative(i.m, k, i.fluid.val[fluid], increment_filter) + if fluid in i.fluid.is_var: + self.jacobian[k, i.fluid.J_col[fluid]] = i.m.val_SI + + k += 1 + + def fluid_dependents(self): + return { + "scalars": [ + [c.m for c in self.inl + self.outl] + for f in self.variable_fluids + ], + "vectors": [{ + c.fluid: set(f) & c.fluid.is_var for c in self.inl + self.outl + } for f in self.variable_fluids] + } + + def energy_balance_func(self): + r""" + Calculate energy balance. + + Returns + ------- + residual : list + Residual value of energy balance. + + .. math:: + + 0 = T_{in} - T_{out,j}\\ + \forall j \in \text{outlets} + """ + residual = [] + T_in = self.inl[0].calc_T() + for o in self.outl: + residual += [T_in - o.calc_T()] + return residual + + def energy_balance_deriv(self, increment_filter, k, dependents=None): + r""" + Calculate partial derivatives of energy balance. + + Parameters + ---------- + increment_filter : ndarray + Matrix for filtering non-changing variables. + + k : int + Position of derivatives in Jacobian matrix (k-th equation). + """ + i = self.inl[0] + dT_dp_in = 0 + dT_dh_in = 0 + if i.p.is_var: + # outlet pressure must be variable as well in this case! + dT_dp_in = dT_mix_dph(i.p.val_SI, i.h.val_SI, i.fluid_data, i.mixing_rule) + if i.h.is_var: + dT_dh_in = dT_mix_pdh(i.p.val_SI, i.h.val_SI, i.fluid_data, i.mixing_rule) + + for o in self.outl: + args = (o.p.val_SI, o.h.val_SI, o.fluid_data, o.mixing_rule) + + dT_dp_out = 0 + if o.p.is_var: + dT_dp_out = -dT_mix_dph(*args) + # pressure is always coupled + self._partial_derivative(i.p, k, dT_dp_in - dT_dp_out) + + if o.h.is_var: + dT_dh_out = -dT_mix_pdh(*args) + + # enthalpy is not necessarily coupled + if i.h._reference_container == o.h._reference_container: + self._partial_derivative(i.h, k, dT_dh_in - dT_dh_out) + else: + self._partial_derivative(i.h, k, dT_dh_in) + self._partial_derivative(o.h, k, dT_dh_out) + + k += 1 + + def energy_balance_dependents(self): + return [ + [self.inl[0].p, self.inl[0].h, o.p, o.h] for o in self.outl + ] diff --git a/src/tespy/connections/__init__.py b/src/tespy/connections/__init__.py index 3741bb556..2eba79104 100644 --- a/src/tespy/connections/__init__.py +++ b/src/tespy/connections/__init__.py @@ -3,4 +3,5 @@ from .bus import Bus # noqa: F401 from .connection import Connection # noqa: F401 from .connection import Ref # noqa: F401 +from .humidairconnection import HAConnection # noqa: F401 from .powerconnection import PowerConnection # noqa: F401 diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index 73d9a9ecc..76eee5f24 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -56,6 +56,7 @@ from tespy.tools.helpers import _is_variable from tespy.tools.helpers import _partial_derivative from tespy.tools.helpers import _partial_derivative_vecvar +from tespy.tools.helpers import seeded_random from tespy.tools.units import SI_UNITS @@ -245,6 +246,9 @@ def _serializable(): def get_variables(self): return {} + def _guess_starting_values(self, units): + pass + def _preprocess(self, row_idx): self.num_eq = 0 @@ -757,7 +761,7 @@ def set_attr(self, **kwargs): # other boolean keywords elif key in ['printout', 'local_design', 'local_offdesign']: if not isinstance(kwargs[key], bool): - msg = ('Please provide the ' + key + ' as boolean.') + msg = f"Please provide the {key} as boolean." logger.error(msg) raise TypeError(msg) else: @@ -887,6 +891,72 @@ def _create_fluid_wrapper(self): fluid, back_end, **wrapper_kwargs ) + def _guess_starting_values(self, units): + # the below part does not work for PowerConnection right now + if sum(self.fluid.val.values()) == 0: + msg = ( + 'The starting value for the fluid composition of the ' + f'connection {self.label} is empty. This might lead to issues ' + 'in the initialisation and solving process as fluid ' + 'property functions can not be called. Make sure you ' + 'specified a fluid composition in all parts of the network.' + ) + logger.warning(msg) + + for key, variable in self.get_variables().items(): + # for connections variables can be presolved and not be var anymore + if variable.is_var: + if not self.good_starting_values: + self._guess_starting_value_from_connected_components(key, units) + + variable.set_SI_from_val0(units) + # variable.set_SI_from_val0() + variable.set_reference_val_SI(variable._val_SI) + + self._precalc_guess_values() + + def _guess_starting_value_from_connected_components(self, key, units): + r""" + Set starting values for fluid properties. + + The component classes provide generic starting values for their inlets + and outlets. + + Parameters + ---------- + c : tespy.connections.connection.Connection + Connection to initialise. + """ + if np.isnan(self.get_attr(key).val0): + # starting value for mass flow is random between 1 and 2 kg/s + # (should be generated based on some hash maybe?) + if key == 'm': + rndm = seeded_random(self.label) + value = float(rndm + 1) + + # generic starting values for pressure and enthalpy + elif key in ['p', 'h']: + # retrieve starting values from component information + val_s = self.source.initialise_source(self, key) + val_t = self.target.initialise_target(self, key) + + if val_s == 0 and val_t == 0: + if key == 'p': + value = 1e5 + elif key == 'h': + value = 1e6 + + elif val_s == 0: + value = val_t + elif val_t == 0: + value = val_s + else: + value = (val_s + val_t) / 2 + + # these values are SI, so they are set to the respective variable + self.get_attr(key).set_reference_val_SI(value) + self.get_attr(key).set_val0_from_SI(units) + def _precalc_guess_values(self): """ Precalculate enthalpy values for connections. @@ -969,6 +1039,7 @@ def _presolve(self): raise TESPyNetworkError(msg) presolved_equations = [] + if self.p.is_set: if self.T_dew.is_set or self.T_bubble.is_set: msg = ( @@ -1240,7 +1311,7 @@ def get_parameters(self): dependents=self.Td_bp_dependents, num_eq=1, quantity="temperature_difference", description="temperature difference to boiling point (deprecated)" - ) + ), } def get_fluid_data(self): @@ -1762,7 +1833,7 @@ def _adjust_to_property_limits(self, nw): self._adjust_enthalpy(fl) # two-phase related - if (self.Td_bp.is_set or self.state.is_set or self.x.is_set or self.td_bubble.is_set or self.td_dew.is_set) and self.it < 10: + if (self.Td_bp.is_set or self.state.is_set or self.x.is_set or self.td_bubble.is_set or self.td_dew.is_set) and self.it < 30: self._adjust_to_two_phase(fl) # mixture diff --git a/src/tespy/connections/humidairconnection.py b/src/tespy/connections/humidairconnection.py new file mode 100644 index 000000000..807196845 --- /dev/null +++ b/src/tespy/connections/humidairconnection.py @@ -0,0 +1,260 @@ +# -*- coding: utf-8 +"""Module of class Connection and class Ref. + + +This file is part of project TESPy (github.com/oemof/tespy). It's copyrighted +by the contributors recorded in the version control history of the file, +available from its original location tespy/connections/humidairconnection.py + +SPDX-License-Identifier: MIT +""" + +import numpy as np +from CoolProp.CoolProp import HAPropsSI + +from tespy.tools import fluid_properties as fp +from tespy.tools.data_containers import FluidComposition as dc_flu +from tespy.tools.data_containers import FluidProperties as dc_prop +from tespy.tools.data_containers import SimpleDataContainer as dc_simple +from tespy.tools.fluid_properties.functions import w_mix_pT_humidair +from tespy.tools.fluid_properties.functions import h_mix_pT +from tespy.tools.fluid_properties.mixtures import _get_fluid_alias, w_mix_fluid_data +from tespy.tools.helpers import seeded_random + +from .connection import Connection +from .connection import connection_registry + + +@connection_registry +class HAConnection(Connection): + + def get_parameters(self): + return { + "fluid": dc_flu( + d=1e-5, + description="mass fractions of the fluid composition" + ), + "m": dc_prop( + quantity="mass_flow", + description="mass flow of dry air (system variable)" + ), + "mHA": dc_prop( + quantity="mass_flow", + description="mass flow of humid air" + ), + "mH2O": dc_prop( + quantity="mass_flow", + description="mass flow of liquid or solid water not contained in humid air" + ), + "p": dc_prop( + quantity="pressure", + description="absolute pressure of the fluid (system variable)" + ), + "h": dc_prop( + quantity="enthalpy", + description="dry air mass specific enthalpy (system variable)" + ), + "w": dc_prop( + quantity="ratio", + description="mass of water per mass of dry air" + ), + "T": dc_prop( + func=self.T_func, + dependents=self.T_dependents, + num_eq=1, + quantity="temperature", + description="temperature of the fluid" + ), + "v": dc_prop( + func=self.v_func, + dependents=self.v_dependents, + num_eq=1, + quantity="volumetric_flow", + description="volumetric flow of the fluid" + ), + "vol": dc_prop( + quantity="specific_volume", + description="specific volume of the fluid (output only)" + ), + "s": dc_prop( + quantity="entropy", + description="specific entropy of the fluid (output only)" + ), + "r": dc_prop( + func=self.r_func, + dependents=self.r_dependents, + num_eq=1, + quantity="ratio", + description="relative humidity" + ), + "fluid_balance": dc_simple( + func=self.fluid_balance_func, + deriv=self.fluid_balance_deriv, + _val=False, num_eq_sets=1, + dependents=self.fluid_balance_dependents, + description="apply an equation which closes the fluid balance with at least two unknown fluid mass fractions" + ), + } + + def _parameter_specification(self, key, value): + if key == "w" or key == "w0": + # specification of w is equivalent to specification of fluid + # composition for humid air + air = 1 / (1 + value) + if key == "w": + self.set_attr(fluid={"air": air, "water": 1 - air}) + else: + self.set_attr(fluid0={"air": air, "water": 1 - air}) + else: + super()._parameter_specification(key, value) + + # for HAConnection mixing rule cannot be modified, is always humidair + def _get_mixing_rule(self): + return "humidair" + + def _set_mixing_rule(self, value): + if value is not None and value != self.mixing_rule: + print(value) + msg = ( + "You cannot change the mixing rule specification for a " + f"Connection of type {self.__class__.__name__}" + ) + raise ValueError(msg) + + mixing_rule = property(_get_mixing_rule, _set_mixing_rule) + + def _guess_starting_values(self, units): + if self.h.is_var and not self.good_starting_values: + value = seeded_random(self.label) + T_rand = 280 + value * (300 - 280) + h = fp.h_mix_pT(1e5, T_rand, self.fluid_data, self.mixing_rule) + self.h.set_reference_val_SI(h) + self._precalc_guess_values() + + def _precalc_guess_values(self): + if not self.h.is_var: + return + + if not self.good_starting_values: + if self.T.is_set: + try: + w = self.calc_w() + self.h.set_reference_val_SI( + HAPropsSI("H", "P", self.p.val_SI, "T", self.T.val_SI, "W", w) + ) + except ValueError: + pass + + def _presolve(self): + + air_alias = _get_fluid_alias("air", self.fluid_data) + water_alias = _get_fluid_alias("water", self.fluid_data) + if not air_alias: + msg = "air must be present in fluid composition" + raise ValueError(msg) + + elif not water_alias: + msg = "water must be present in fluid composition" + raise ValueError(msg) + + if len(self.fluid.is_var) > 0: + return [] + + presolved_equations = [] + + if self.h.is_var and not self.p.is_var: + if self.T.is_set: + self.h.set_reference_val_SI(h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule)) + self.h._potential_var = False + if "T" in self._equation_set_lookup.values(): + presolved_equations += ["T"] + + presolved_equations = [ + key for parameter in presolved_equations + for key, value in self._equation_set_lookup.items() + if value == parameter + ] + return presolved_equations + + def _adjust_to_property_limits(self, nw): + + if self.p.is_var: + if self.p.val_SI < 100: + self.p.val_SI = 101 + elif self.p.val_SI > 100e5: + self.p.val_SI = 99e5 + + if self.h.is_var: + # TODO: check minimum temperature how it matches minimum humidity ratio + d = self.h._reference_container._d + hmin = HAPropsSI("H", "T", -50 + 273.15, "P", self.p.val_SI, "R", 1) + if self.h.val_SI < hmin: + delta = max(abs(self.h.val_SI * d), d) * 5 + self.h.set_reference_val_SI(hmin + delta) + + else: + # TODO: where to get reasonable hmax from?! + hmax = HAPropsSI("H", "T", 300 + 273.15, "P", self.p.val_SI, "R", 0) + if self.h.val_SI > hmax: + delta = max(abs(self.h.val_SI * d), d) * 5 + self.h.set_reference_val_SI(hmax - delta) + + @classmethod + def _result_attributes(cls): + return ["m", "mHA", "mH2O", "p", "h", "T", "w", "s", "vol", "v", "r"] + + @classmethod + def _print_attributes(cls): + return ["m", "mHA", "mH2O", "p", "h", "T", "w", "r"] + + def calc_r(self): + w = self.calc_w() + try: + return HAPropsSI("R", "P", self.p.val_SI, "T", self.T.val_SI, "W", w) + except ValueError as e: + value = str(e).split("value (")[1].split(")")[0] + return float(value) + + def r_func(self): + return self.r.val_SI - self.calc_r() + + def r_dependents(self): + water_alias = _get_fluid_alias("H2O", self.fluid_data) + # water alias is already a set + return { + "scalars": [self.p, self.h], + "vectors": [{self.fluid: water_alias}] + } + + def calc_w(self): + return w_mix_pT_humidair(self.p.val_SI, self.T.val_SI, self.fluid_data) + + def calc_results(self, units): + self.T.val_SI = self.calc_T() + self.vol.val_SI = self.calc_vol() # Mixture volume per mass of dry air + self.v.val_SI = self.vol.val_SI * self.m.val_SI # Mixture volume flow rate + # handle the water fraction + self.w.val_SI = self.calc_w() + # # Convert from kg water/kg dry air to mass fraction of water in humid air + # x_h2o = self.w.val_SI / (1 + self.w.val_SI) + # x_air = 1 - x_h2o + self.mHA.val_SI = self.m.val_SI * (1 + self.w.val_SI) + w_mixture = w_mix_fluid_data(self.fluid_data) + # Calculate kg water/kg dry air that is not in the humid air + delta_w = w_mixture - self.w.val_SI + self.mH2O.val_SI = self.m.val_SI * delta_w + self.r.val_SI = self.calc_r() + # if self.r.val_SI > 1: + # self.r.val_SI = np.nan + + for prop in self._result_attributes(): + param = self.get_attr(prop) + if not param.is_set: + param.set_val_from_SI(units) + + self.m.set_val0_from_SI(units) + self.p.set_val0_from_SI(units) + self.h.set_val0_from_SI(units) + self.fluid.val0 = self.fluid.val.copy() + + return True diff --git a/src/tespy/connections/powerconnection.py b/src/tespy/connections/powerconnection.py index 3503e346a..7589c0a41 100644 --- a/src/tespy/connections/powerconnection.py +++ b/src/tespy/connections/powerconnection.py @@ -130,6 +130,10 @@ def set_attr(self, **kwargs): logger.error(msg) raise KeyError(msg) + def _guess_starting_values(self, units): + if self.E.is_var and not self.good_starting_values: + self.E.set_reference_val_SI(0.0) + def _precalc_guess_values(self): pass diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index 39cd92667..d2e574adc 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -16,6 +16,7 @@ import json import math import os +from pathlib import Path import warnings from time import time @@ -1764,6 +1765,26 @@ def _prepare_design(self): data = _local_designs[path][c] # write data self._write_design_state_to_component(cp, data) + # this is a hack to write the connection specifications of the + # respective component into the component. In the component + # itself, the value from here should be utilized instead of the + # .design value of the connection in case local_offdesign is + # set to True + if cp.inl[0].label in _local_designs[path][cp.inl[0].__class__.__name__].index: + cp._connection_offdesign = { + c.label: _local_designs[path][c.__class__.__name__].loc[c.label] + for c in cp.inl + cp.outl + } + else: + # remap to actual inlet and outlet labels + cp._connection_offdesign = { + c.label: _local_designs[path][c.__class__.__name__].loc[f"hx_in{k + 1}"] + for k, c in enumerate(cp.inl) + } + cp._connection_offdesign.update({ + c.label: _local_designs[path][c.__class__.__name__].loc[f"hx_out{k + 1}"] + for k, c in enumerate(cp.outl) + }) # unset design parameters for var in cp.design: @@ -2030,27 +2051,7 @@ def _set_starting_values(self): df = dfs[c.__class__.__name__] self._write_starting_values_to_connection(c, df) - if type(c) == Connection: - # the below part does not work for PowerConnection right now - if sum(c.fluid.val.values()) == 0: - msg = ( - 'The starting value for the fluid composition of the ' - f'connection {c.label} is empty. This might lead to issues ' - 'in the initialisation and solving process as fluid ' - 'property functions can not be called. Make sure you ' - 'specified a fluid composition in all parts of the network.' - ) - logger.warning(msg) - - for key, variable in c.get_variables().items(): - # for connections variables can be presolved and not be var anymore - if variable.is_var: - if not c.good_starting_values: - self._guess_starting_value_from_connected_components(c, key) - - variable.set_SI_from_val0(self.units) - # variable.set_SI_from_val0() - variable.set_reference_val_SI(variable._val_SI) + c._guess_starting_values(self.units) for cp in self.comps["object"]: for key, variable in cp.get_variables().items(): @@ -2061,70 +2062,35 @@ def _set_starting_values(self): variable.set_SI_from_val(self.units) variable.set_reference_val_SI(variable._val_SI) - for c in self.conns['object']: - c._precalc_guess_values() - msg = 'Generic fluid property specification complete.' logger.debug(msg) - def _guess_starting_value_from_connected_components(self, c, key): - r""" - Set starting values for fluid properties. - - The component classes provide generic starting values for their inlets - and outlets. - - Parameters - ---------- - c : tespy.connections.connection.Connection - Connection to initialise. - """ - if np.isnan(c.get_attr(key).val0): - # starting value for mass flow is random between 1 and 2 kg/s - # (should be generated based on some hash maybe?) - if key == 'm': - seed = abs(hash(c.label)) % (2**32) - rng = np.random.default_rng(seed=seed) - value = float(rng.random() + 1) - - # generic starting values for pressure and enthalpy - elif key in ['p', 'h']: - # retrieve starting values from component information - val_s = c.source.initialise_source(c, key) - val_t = c.target.initialise_target(c, key) - - if val_s == 0 and val_t == 0: - if key == 'p': - value = 1e5 - elif key == 'h': - value = 1e6 - - elif val_s == 0: - value = val_t - elif val_t == 0: - value = val_s - else: - value = (val_s + val_t) / 2 - - elif key == 'E': - value = 0.0 - - # these values are SI, so they are set to the respective variable - c.get_attr(key).set_reference_val_SI(value) - c.get_attr(key).set_val0_from_SI(self.units) @staticmethod - def _load_network_state(json_path): + def _load_network_state(json_path: str | bytes | bytearray | Path): r""" Read network state from given file. Parameters ---------- - json_path : str + json_path : str | bytes | bytearray | Path Path to network information. """ - with open(json_path, "r") as f: - data = json.load(f) + data = None + if not isinstance(json_path, Path): + try: + data = json.loads(json_path) + except json.JSONDecodeError as e: + msg = ( + "The provided json_path could not be decoded. If this is not " + "a valid json string, please provide a valid file path instead of " + "%s" + ) + logger.debug(msg, str(json_path)) + pass + if data is None: + with open(json_path, "r") as f: + data = json.load(f) dfs = {} if "Connection" in data["Connection"]: @@ -2365,6 +2331,23 @@ def _get_linear_dependents_by_variable_index(self, idx) -> list: variables = [self._variable_lookup[v] for v in dependents["variables"]] variable_list = [(v["object"].label, v["property"]) for v in variables] return variable_list + + def get_sorted_residual_index(self) -> list[int]: + """Get the sorted array of residual indices. + + Returns + ------- + list[int] + List of variable numbers, the index values. + """ + # vars: dict[tuple[int, str], dict] = self.get_variables() + sidx: list[int] = list(np.argsort(np.abs(self.residual))[::-1]) + # sres = np.array([self.residual[i] for i in sidx]) + # chis = self.residual_history.shape[1] + # for i in range(2, n): + # sres = np.vstack((sres, [self.residual_history[i-2][j] for j in sidx])) + # sres = np.vstack((sres, self.residual_history[-n+1:, :][:, sidx].T)) + return sidx def solve(self, mode, init_path=None, design_path=None, max_iter=50, min_iter=4, init_only=False, init_previous=True, @@ -2852,18 +2835,19 @@ def _update_variables(self): # get_J_col yet relax = 1 if self.robust_relax: - if self.iter < 3: - relax = 0.25 - elif self.iter < 5: - relax = 0.5 - elif self.iter < 8: - relax = 0.75 + # relax_values = [ + # (0.05, 0.05 * self.max_iter), + # (0.10, 0.10 * self.max_iter), + # (0.25, 0.25 * self.max_iter), + # (0.50, 0.50 * self.max_iter) + # ] + relax = 0.05 + 0.95 * min(1, self.iter / (0.25 * self.max_iter)) for _, data in self.variables_dict.items(): if data["variable"] in ["m", "h", "E"]: container = data["obj"] container._val_SI += increment[container.J_col] * relax - elif data["variable"] == "p": + elif data["variable"] in ["p", "w"]: container = data["obj"] p_relax = max( 1, -2 * increment[container.J_col] / container.val_SI @@ -3556,18 +3540,26 @@ def export(self, json_file_path=None): return export - def save(self, json_file_path): + def save(self, json_file_path: str | Path | None) -> None | str: r""" Dump the results to a json style output. Parameters ---------- - json_file_path : str + json_file_path : str | Path | None Filename to dump results into. + Returns + ------- + None + If a file path is provided, results are saved to file. + str + If no file path is provided, results are returned as string. + Note ---- - Results will be saved to specified file path + Results will be saved to specified file path in json format. If no + file path is provided, the results will be returned as string. """ dump = {} @@ -3578,6 +3570,9 @@ def save(self, json_file_path): dump = hlp._nested_dict_of_dataframes_to_dict(dump) + if json_file_path is None: + return json.dumps(dump, indent=2) + with open(json_file_path, "w") as f: json.dump(dump, f) diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py index 00a0bf3c4..7b19da5bd 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/src/tespy/tools/fluid_properties/functions.py @@ -10,6 +10,7 @@ SPDX-License-Identifier: MIT """ +from CoolProp.CoolProp import HAPropsSI from tespy.tools.global_vars import FLUID_ALIASES from tespy.tools.logger import logger @@ -25,6 +26,7 @@ from .mixtures import T_MIX_PS_REVERSE from .mixtures import V_MIX_PT_DIRECT from .mixtures import VISCOSITY_MIX_PT_DIRECT +from .mixtures import w_mix_pT_humidair, w_mix_ph_humidair, w_mix_ps_humidair def isentropic(p_1, h_1, p_2, fluid_data, mixing_rule=None, T0=None): @@ -138,11 +140,15 @@ def T_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): return pure_fluid["wrapper"].T_ph(p, h) else: _check_mixing_rule(mixing_rule, T_MIX_PH_REVERSE, "temperature (from enthalpy)") - kwargs = { - "p": p, "target_value": h, "fluid_data": fluid_data, "T0": T0, - "f": T_MIX_PH_REVERSE[mixing_rule] - } - return inverse_temperature_mixture(**kwargs) + if mixing_rule == "humidair": + w = w_mix_ph_humidair(p, h, fluid_data) + return HAPropsSI("T", "P", p, "H", h, "W", w) + else: + kwargs = { + "p": p, "target_value": h, "fluid_data": fluid_data, "T0": T0, + "f": T_MIX_PH_REVERSE[mixing_rule] + } + return inverse_temperature_mixture(**kwargs) def dT_mix_pdh(p, h, fluid_data, mixing_rule=None, T0=None): @@ -296,12 +302,16 @@ def T_mix_ps(p, s, fluid_data, mixing_rule=None, T0=None): pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].T_ps(p, s) else: - _check_mixing_rule(mixing_rule, T_MIX_PS_REVERSE, "temperature (from entropy)") - kwargs = { - "p": p, "target_value": s, "fluid_data": fluid_data, "T0": T0, - "f": T_MIX_PS_REVERSE[mixing_rule] - } - return inverse_temperature_mixture(**kwargs) + if mixing_rule == "humidair": + w = w_mix_ps_humidair(p, s, fluid_data) + return HAPropsSI("T", "P", p, "S", s, "W", w) + else: + _check_mixing_rule(mixing_rule, T_MIX_PS_REVERSE, "temperature (from entropy)") + kwargs = { + "p": p, "target_value": s, "fluid_data": fluid_data, "T0": T0, + "f": T_MIX_PS_REVERSE[mixing_rule] + } + return inverse_temperature_mixture(**kwargs) def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): @@ -309,8 +319,12 @@ def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): pure_fluid = get_pure_fluid(fluid_data) return 1 / pure_fluid["wrapper"].d_ph(p, h) else: - T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) - return v_mix_pT(p, T, fluid_data, mixing_rule) + if mixing_rule == "humidair": + w = w_mix_ph_humidair(p, h, fluid_data) + return HAPropsSI("V", "P", p, "H", h, "W", w) + else: + T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) + return v_mix_pT(p, T, fluid_data, mixing_rule) def dv_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None): @@ -341,8 +355,12 @@ def viscosity_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].viscosity_ph(p, h) else: - T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) - return viscosity_mix_pT(p, T, fluid_data, mixing_rule) + if mixing_rule == "humidair": + w = w_mix_ph_humidair(p, h, fluid_data) + return HAPropsSI("Visc", "P", p, "H", h, "W", w) + else: + T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) + return viscosity_mix_pT(p, T, fluid_data, mixing_rule) def viscosity_mix_pT(p, T, fluid_data, mixing_rule=None): diff --git a/src/tespy/tools/fluid_properties/helpers.py b/src/tespy/tools/fluid_properties/helpers.py index 77dc3700e..7b2f7f8f1 100644 --- a/src/tespy/tools/fluid_properties/helpers.py +++ b/src/tespy/tools/fluid_properties/helpers.py @@ -49,6 +49,9 @@ def get_pure_fluid(fluid_data): def single_fluid(fluid_data): r"""Return the name of the pure fluid in a fluid vector.""" + if "_HUMID_AIR" in fluid_data: + return None + if get_number_of_fluids(fluid_data) > 1: return None else: diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index d26cdcb10..8b1a42139 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -13,6 +13,8 @@ import math +from CoolProp.CoolProp import HAPropsSI + from tespy.tools.global_vars import FLUID_ALIASES from tespy.tools.global_vars import gas_constants @@ -36,7 +38,7 @@ def h_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): def h_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) mass_fractions_gas, molar_fraction_gas, mass_liquid, _, p_sat, pp_water = cond_check(p, T, fluid_data, water_alias) @@ -88,6 +90,56 @@ def h_mix_pT_incompressible(p, T, fluid_data, **kwargs): return h +def w_mix_fluid_data(fluid_data): + + water_alias = _get_fluid_alias("H2O", fluid_data) + water_alias = next(iter(water_alias)) + + air_alias = _get_fluid_alias("air", fluid_data) + air_alias = next(iter(air_alias)) + + return ( + fluid_data[water_alias]["mass_fraction"] + / fluid_data[air_alias]["mass_fraction"] + ) + +def w_mix_pTrh_humidair(p, T, rh): + return HAPropsSI("W", "P", p, "T", T, "RH", rh) # kg water/kg dry air + +def w_mix_pT_humidair(p, T, fluid_data, **kwargs): + w_def = w_mix_fluid_data(fluid_data) + w_max = w_mix_pTrh_humidair(p, T, 1.0) + if w_def > w_max: + _msg = f"Humidity ratio {w_def:.4f} exceeds maximum value of {w_max:.4f} for given p and T. Check fluid composition." + return w_max + return w_def + +def h_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_pT_humidair(p, T, fluid_data, **kwargs) + return HAPropsSI("H", "P", p, "T", T, "W", w) + +def w_mix_phrh_humidair(p, h, rh): + return HAPropsSI("W", "P", p, "H", h, "RH", rh) # kg water/kg dry air + +def w_mix_ph_humidair(p, h, fluid_data, **kwargs): + w_def = w_mix_fluid_data(fluid_data) + w_max = w_mix_phrh_humidair(p, h, 1.0) + if w_def > w_max: + _msg = f"Humidity ratio {w_def:.4f} exceeds maximum value of {w_max:.4f} for given p and T. Check fluid composition." + return w_max + return w_def + +def w_mix_psrh_humidair(p, s, rh): + return HAPropsSI("W", "P", p, "S", s, "RH", rh) # kg water/kg dry air + +def w_mix_ps_humidair(p, s, fluid_data, **kwargs): + w_def = w_mix_fluid_data(fluid_data) + w_max = w_mix_psrh_humidair(p, s, 1.0) + if w_def > w_max: + _msg = f"Humidity ratio {w_def:.4f} exceeds maximum value of {w_max:.4f} for given p and T. Check fluid composition." + return w_max + return w_def + def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): molar_fractions = get_molar_fractions(fluid_data) @@ -103,7 +155,7 @@ def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): def s_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) mass_fractions_gas, molar_fraction_gas, mass_liquid, _, p_sat, pp_water = cond_check(p, T, fluid_data, water_alias) @@ -136,6 +188,11 @@ def s_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return s +def s_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_fluid_data(fluid_data) + return HAPropsSI("S", "P", p, "T", T, "W", w) + + def v_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): molar_fractions = get_molar_fractions(fluid_data) @@ -151,7 +208,7 @@ def v_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): def v_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) _, molar_fraction_gas, mass_liquid, _, p_sat, pp_water = cond_check(p, T, fluid_data, water_alias) @@ -183,6 +240,11 @@ def v_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return v +def v_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_fluid_data(fluid_data) + return HAPropsSI("V", "P", p, "T", T, "W", w) + + def viscosity_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): r""" Calculate dynamic viscosity from pressure and temperature. @@ -239,10 +301,15 @@ def viscosity_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return viscosity +def viscosity_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_fluid_data(fluid_data) + return HAPropsSI("Visc", "P", p, "T", T, "W", w) + + def exergy_chemical_ideal_cond(pamb, Tamb, fluid_data, Chem_Ex): molar_fractions = get_molar_fractions(fluid_data) - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) _, molar_fractions_gas, _, molar_liquid, _, _ = cond_check( @@ -277,9 +344,9 @@ def exergy_chemical_ideal_cond(pamb, Tamb, fluid_data, Chem_Ex): return ex_chemical * 1e3 # Data from Chem_Ex are in kJ / mol -def _water_in_mixture(fluid_data): +def _get_fluid_alias(fluid, fluid_data): return ( - FLUID_ALIASES.get_fluid("H2O") + FLUID_ALIASES.get_fluid(fluid) & set([ f for f in fluid_data if _is_larger_than_precision(fluid_data[f]["mass_fraction"]) @@ -351,14 +418,16 @@ def cond_check(p, T, fluid_data, water_alias): T_MIX_PH_REVERSE = { "ideal": h_mix_pT_ideal, "ideal-cond": h_mix_pT_ideal_cond, - "incompressible": h_mix_pT_incompressible + "incompressible": h_mix_pT_incompressible, + "humidair": h_mix_pT_humidair } T_MIX_PS_REVERSE = { "ideal": s_mix_pT_ideal, "ideal-cond": s_mix_pT_ideal_cond, - "incompressible": s_mix_pT_incompressible + "incompressible": s_mix_pT_incompressible, + "humidair": s_mix_pT_humidair } @@ -366,28 +435,32 @@ def cond_check(p, T, fluid_data, water_alias): "ideal": h_mix_pT_ideal, "ideal-cond": h_mix_pT_ideal_cond, "incompressible": h_mix_pT_incompressible, - "forced-gas": h_mix_pT_forced_gas + "forced-gas": h_mix_pT_forced_gas, + "humidair": h_mix_pT_humidair } S_MIX_PT_DIRECT = { "ideal": s_mix_pT_ideal, "ideal-cond": s_mix_pT_ideal_cond, - "incompressible": s_mix_pT_incompressible + "incompressible": s_mix_pT_incompressible, + "humidair": s_mix_pT_humidair } V_MIX_PT_DIRECT = { "ideal": v_mix_pT_ideal, "ideal-cond": v_mix_pT_ideal_cond, - "incompressible": v_mix_pT_incompressible + "incompressible": v_mix_pT_incompressible, + "humidair": v_mix_pT_humidair } VISCOSITY_MIX_PT_DIRECT = { "ideal": viscosity_mix_pT_ideal, "ideal-cond": viscosity_mix_pT_ideal, - "incompressible": viscosity_mix_pT_incompressible + "incompressible": viscosity_mix_pT_incompressible, + "humidair": viscosity_mix_pT_humidair } EXERGY_CHEMICAL = { diff --git a/src/tespy/tools/helpers.py b/src/tespy/tools/helpers.py index 92cc76478..1a114d8e5 100644 --- a/src/tespy/tools/helpers.py +++ b/src/tespy/tools/helpers.py @@ -16,6 +16,7 @@ from copy import deepcopy import pandas as pd +import numpy as np from tespy import __datapath__ from tespy.tools import logger @@ -700,6 +701,23 @@ def central_difference(function=None, parameter=None, delta=None, **kwargs): return (function(**upper) - function(**lower)) / (2 * delta) +def seeded_random_generator(seed_value: str | bytes | int) -> np.random.Generator: + """Generate a reproducible random number generator based on a seed value.""" + if isinstance(seed_value, str): + seed_value = seed_value.encode("utf-8") + if isinstance(seed_value, bytes): + seed_value = int.from_bytes(seed_value, "little", signed=False) % (2**32) + if not isinstance(seed_value, int): + raise ValueError("Seed value must be of type str, bytes or int.") + return np.random.default_rng(seed=seed_value) + + +def seeded_random(seed_value: str | bytes | int) -> float: + """Generate a reproducible random number between 0 and 1 based on a seed value.""" + rng = seeded_random_generator(seed_value) + return rng.random() + + def get_basic_path(): """ Return the basic tespy path and creates it if necessary. diff --git a/tests/test_connections.py b/tests/test_connections.py index 86401678a..1a3702047 100644 --- a/tests/test_connections.py +++ b/tests/test_connections.py @@ -29,6 +29,7 @@ from tespy.connections import Ref from tespy.connections.connection import ConnectionBase from tespy.connections.connection import connection_registry +from tespy.connections.humidairconnection import HAConnection from tespy.networks import Network from tespy.tools.data_containers import FluidProperties as dc_prop from tespy.tools.fluid_properties.functions import T_bubble_p @@ -516,7 +517,7 @@ def test_all_classes_in_registry(obj): def make_connection(cls): - if cls == Connection: + if cls == Connection or cls == HAConnection: return cls(Source(""), "out1", Sink(""), "in1") elif cls == PowerConnection: return cls(PowerSource(""), "power", PowerSink(""), "power") diff --git a/tests/test_tools/test_fluid_properties/humidair.ipynb b/tests/test_tools/test_fluid_properties/humidair.ipynb new file mode 100644 index 000000000..d93eef43a --- /dev/null +++ b/tests/test_tools/test_fluid_properties/humidair.ipynb @@ -0,0 +1,187 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "5796d5ee", + "metadata": {}, + "outputs": [], + "source": [ + "from tespy.connections import HAConnection, Connection, Ref\n", + "from tespy.components import Source, Sink, MovingBoundaryHeatExchanger\n", + "from tespy.networks import Network\n", + "\n", + "from CoolProp.CoolProp import HAPropsSI" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "947c79c2", + "metadata": {}, + "outputs": [], + "source": [ + "def get_water_air_mixture_fractions_pTR(p, T, R):\n", + " w = HAPropsSI(\"W\", \"P\", p, \"T\", T, \"R\", R)\n", + " return get_water_air_mixture_fractions_w(w)\n", + "\n", + "\n", + "def get_water_air_mixture_fractions_w(w):\n", + " air = 1 / (1 + w)\n", + " return {\"Air\": air, \"water\": 1 - air}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7351d18", + "metadata": {}, + "outputs": [], + "source": [ + "from tespy.tools.fluid_properties.functions import _get_humid_air_humidity_ratio\n", + "\n", + "\n", + "def _calculate_dQ_freezing_water(T_evap, T_air, air_fluid_data, air_pressure, air_mass_flow):\n", + " dh_water = 4000 * (T_air - 273.15)\n", + " dh_freeze = 333000\n", + " dh_ice = 2000 * (273.15 - T_evap)\n", + " w_actual = _get_humid_air_humidity_ratio(air_fluid_data)\n", + " w_saturated = HAPropsSI(\"W\", \"P\", air_pressure, \"T\", T_air, \"R\", 1)\n", + "\n", + " m_water_freeze = (w_actual - w_saturated) * air_mass_flow\n", + " return (dh_water + dh_freeze + dh_ice) * m_water_freeze\n", + "\n", + "\n", + "class MovingBoundaryHeatExchangerWithFreeze(MovingBoundaryHeatExchanger):\n", + "\n", + " def energy_balance_func(self):\n", + " residual = super().energy_balance_func()\n", + " T_evap = self.outl[1].calc_T()\n", + " if T_evap < 273.15:\n", + " T_air_out = self.outl[0].calc_T()\n", + " if T_air_out > 273.15:\n", + " dQ_freeze = _calculate_dQ_freezing_water(\n", + " T_evap, T_air_out, self.outl[0].fluid_data, self.outl[0].p.val_SI, self.inl[0].m.val_SI\n", + " )\n", + " residual -= dQ_freeze\n", + "\n", + " return residual\n", + "\n", + " def energy_balance_dependents(self):\n", + " return [\n", + " self.inl[0].m, self.inl[0].h,\n", + " self.outl[0].p, self.outl[0].h,\n", + " self.inl[1].m, self.inl[1].h,\n", + " self.outl[1].p, self.outl[1].p\n", + " ]\n", + "\n", + " def calc_parameters(self):\n", + " super().calc_parameters()\n", + " self.Q.val_SI = -self.inl[1].m.val_SI * (self.outl[1].h.val_SI - self.inl[1].h.val_SI)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10a7b936", + "metadata": {}, + "outputs": [], + "source": [ + "nw = Network()\n", + "nw.units.set_defaults(\n", + " temperature=\"°C\",\n", + " pressure=\"bar\",\n", + " heat=\"kW\"\n", + ")\n", + "\n", + "\n", + "so = Source(\"source\")\n", + "hex = MovingBoundaryHeatExchanger(\"heat exchanger\")\n", + "si = Sink(\"sink\")\n", + "\n", + "so2 = Source(\"refrigerant source\")\n", + "si2 = Sink(\"refrigerant sink\")\n", + "\n", + "a1 = Connection(so2, \"out1\", hex, \"in2\", label=\"a1\")\n", + "a2 = Connection(hex, \"out2\", si2, \"in1\", label=\"a2\")\n", + "\n", + "c1 = HAConnection(so, \"out1\", hex, \"in1\", label=\"c1\")\n", + "c2 = HAConnection(hex, \"out1\", si, \"in1\", label=\"c2\")\n", + "\n", + "nw.add_conns(a1, a2, c1, c2)\n", + "\n", + "a1.set_attr(fluid={\"NH3\": 1}, m=2, x=0.3)\n", + "a2.set_attr(td_dew=5, T_dew=-20)\n", + "\n", + "c1.set_attr(\n", + " p=1,\n", + " T=15,\n", + " w0=0.004, # -> implicitly sets starting value for composition\n", + " #w=0.004, -> use this to directly impose fluid composition, then do not specify r and fluid balance for this example\n", + " r=1,\n", + " fluid_balance=True,\n", + ")\n", + "c2.set_attr(T=5)\n", + "\n", + "hex.set_attr(dp1=0, dp2=0)\n", + "\n", + "nw.solve(\"design\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3fc42bdf", + "metadata": {}, + "outputs": [], + "source": [ + "nw.print_results()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd871dc0", + "metadata": {}, + "outputs": [], + "source": [ + "heat, T_hot, T_cold, _, _ = hex.calc_sections()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f9ed087", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "\n", + "\n", + "plt.plot(heat, T_hot)\n", + "plt.plot(heat, T_cold)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "tespy", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}