From e57e17981e0ff182026db2745a382638d09963b5 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Mon, 16 Oct 2023 23:32:30 +0200 Subject: [PATCH 01/22] Add idea for the implementation of solutions --- src/tespy/connections/connection.py | 87 +++++++++++--------- src/tespy/tools/fluid_properties/helpers.py | 2 +- src/tespy/tools/fluid_properties/mixtures.py | 48 ++++++++++- src/tespy/tools/fluid_properties/wrappers.py | 7 +- 4 files changed, 100 insertions(+), 44 deletions(-) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index ceb64a66f..5089f06bf 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -719,6 +719,8 @@ def build_fluid_data(self): "mass_fraction": self.fluid.val[fluid] } for fluid in self.fluid.val } + if self.mixing_rule == "incomp-solution": + self.fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([self.fluid.val["LiBr"]]) def primary_ref_func(self, k, **kwargs): variable = kwargs["variable"] @@ -902,50 +904,53 @@ def solve(self, increment_filter): data.deriv(k, **data.func_params) def calc_results(self): - self.T.val_SI = self.calc_T() - number_fluids = get_number_of_fluids(self.fluid_data) _converged = True - if number_fluids > 1: - h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule) - if abs(h_from_T - self.h.val_SI) > ERR ** .5: - self.T.val_SI = np.nan - self.vol.val_SI = np.nan - self.v.val_SI = np.nan - self.s.val_SI = np.nan - msg = ( - "Could not find a feasible value for mixture temperature at " - f"connection {self.label}. The values for temperature, " - "specific volume, volumetric flow and entropy are set to nan." - ) - logger.error(msg) - _converged = False - - else: - try: - if not self.x.is_set: - self.x.val_SI = self.calc_x() - except ValueError: - self.x.val_SI = np.nan - try: - if not self.Td_bp.is_set: - self.Td_bp.val_SI = self.calc_Td_bp() - except ValueError: - self.x.val_SI = np.nan - - if _converged: - self.vol.val_SI = self.calc_vol() - self.v.val_SI = self.vol.val_SI * self.m.val_SI - self.s.val_SI = self.calc_s() + try: + self.T.val_SI = self.calc_T() + number_fluids = get_number_of_fluids(self.fluid_data) + if number_fluids > 1: + h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule) + if abs(h_from_T - self.h.val_SI) > ERR ** .5: + self.T.val_SI = np.nan + self.vol.val_SI = np.nan + self.v.val_SI = np.nan + self.s.val_SI = np.nan + msg = ( + "Could not find a feasible value for mixture temperature at " + f"connection {self.label}. The values for temperature, " + "specific volume, volumetric flow and entropy are set to nan." + ) + logger.error(msg) + _converged = False - for prop in fpd.keys(): - self.get_attr(prop).val = convert_from_SI( - prop, self.get_attr(prop).val_SI, self.get_attr(prop).unit - ) + else: + try: + if not self.x.is_set: + self.x.val_SI = self.calc_x() + except ValueError: + self.x.val_SI = np.nan + try: + if not self.Td_bp.is_set: + self.Td_bp.val_SI = self.calc_Td_bp() + except ValueError: + self.x.val_SI = np.nan + + if _converged: + self.vol.val_SI = self.calc_vol() + self.v.val_SI = self.vol.val_SI * self.m.val_SI + self.s.val_SI = self.calc_s() + + for prop in fpd.keys(): + self.get_attr(prop).val = convert_from_SI( + prop, self.get_attr(prop).val_SI, self.get_attr(prop).unit + ) - self.m.val0 = self.m.val - self.p.val0 = self.p.val - self.h.val0 = self.h.val - self.fluid.val0 = self.fluid.val.copy() + self.m.val0 = self.m.val + self.p.val0 = self.p.val + self.h.val0 = self.h.val + self.fluid.val0 = self.fluid.val.copy() + except (ValueError, KeyError): + pass def check_pressure_bounds(self, fluid): if self.p.val_SI > self.fluid.wrapper[fluid]._p_max: diff --git a/src/tespy/tools/fluid_properties/helpers.py b/src/tespy/tools/fluid_properties/helpers.py index 27740208a..226648579 100644 --- a/src/tespy/tools/fluid_properties/helpers.py +++ b/src/tespy/tools/fluid_properties/helpers.py @@ -69,7 +69,7 @@ def inverse_temperature_mixture(p=None, target_value=None, fluid_data=None, T0=N valmin, valmax = get_mixture_temperature_range(fluid_data) if T0 is None: T0 = (valmin + valmax) / 2.0 - + T0 = 320 function_kwargs = { "p": p, "fluid_data": fluid_data, "T": T0, "function": f, "parameter": "T" , "delta": 0.01 diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index b156a133f..0cdc4c240 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -34,6 +34,38 @@ def h_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): return h +def xsat_pT_incomp_solution(p=None, T=None, fluid_data=None, **kwargs): + X_min = 0 + X_max = 0.75 + X = .3 + d = 1e-5 + iter = 0 + while True: + res = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X) - p + upper = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X + d) + lower = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X - d) + + deriv = (upper - lower) / (2 * d) + X -= res / deriv + + if abs(res) < 1e-6: + if X < X_min: + return X_min + elif X > X_max: + return X_max + break + + if X >= X_max: + X = X_max - d + elif X <= X_min: + X = X_min + d + + iter += 1 + if iter > 10: + break + return X + + def h_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): water_alias = _water_in_mixture(fluid_data) @@ -85,6 +117,19 @@ def h_mix_pT_incompressible(p, T, fluid_data, **kwargs): return h +def h_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): + if p < fluid_data["LiBr"]["wrapper"].p_sat(T): + x_old = fluid_data["LiBr"]["mass_fraction"] + x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs) + x_water_gas = x_new - x_old + h = fluid_data["LiBr"]["wrapper"].h_pT(p, T - .005) * (1 - x_water_gas) + h += fluid_data["water"]["wrapper"].h_QT(1, T) * x_water_gas + fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([x_old]) + return h + else: + return fluid_data["LiBr"]["wrapper"].h_pT(p, T) + + def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): molar_fractions = get_molar_fractions(fluid_data) @@ -336,7 +381,8 @@ 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, + "incomp-solution": h_mix_pT_incompressible_solution } diff --git a/src/tespy/tools/fluid_properties/wrappers.py b/src/tespy/tools/fluid_properties/wrappers.py index 97ce00647..01099033e 100644 --- a/src/tespy/tools/fluid_properties/wrappers.py +++ b/src/tespy/tools/fluid_properties/wrappers.py @@ -133,7 +133,7 @@ def _set_constants(self): self._p_min = 1e2 self._p_max = 1e8 self._p_crit = 1e8 - self._T_crit = None + self._T_crit = self._T_max self._molar_mass = 1 try: # how to know that we have a binary mixture? @@ -164,6 +164,11 @@ def get_T_max(self, p): def isentropic(self, p_1, h_1, p_2): return self.h_ps(p_2, self.s_ph(p_1, h_1)) + def psat_Tx(self, T, x): + self.AS.set_mass_fractions([x]) + self.AS.update(CP.QT_INPUTS, 0, T) + return self.AS.p() + def T_ph(self, p, h): self.AS.update(CP.HmassP_INPUTS, h, p) return self.AS.T() From 385b92ef3691c22537425ace1d6a4d478e720638 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Mon, 16 Oct 2023 23:32:48 +0200 Subject: [PATCH 02/22] Add draft for absorber component --- test_libr.py | 335 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 335 insertions(+) create mode 100644 test_libr.py diff --git a/test_libr.py b/test_libr.py new file mode 100644 index 000000000..f10a64c20 --- /dev/null +++ b/test_libr.py @@ -0,0 +1,335 @@ +from CoolProp.CoolProp import AbstractState +import CoolProp as CP +import numpy as np + +from tespy.components.component import Component + + +libr = AbstractState("INCOMP", "LiBr") + +libr.set_mass_fractions([0.5]) +libr.update(CP.QT_INPUTS, 0, 300) +libr.p() +libr.set_mass_fractions([0.2]) +libr.update(CP.QT_INPUTS, 0, 300) +libr.p() + +# saturation pressure by temperature and libr mass fraction x + +def psat_TX(libr, T, X): + libr.set_mass_fractions([X]) + libr.update(CP.QT_INPUTS, 0, T) + return libr.p() + + +def Xsat_Tp(libr, T, p): + X_min = 0 + X_max = 0.75 + X = .3 + d = 1e-5 + while True: + res = psat_TX(libr, T, X) - p + deriv = (psat_TX(libr, T, X + d) - psat_TX(libr, T, X - d)) / (2 * d) + X -= res / deriv + # print(X) + + if abs(res) < 1e-6: + if X < X_min: + return X_min + elif X > X_max: + return X_max + break + + if X >= X_max: + X = X_max - d + elif X <= X_min: + X = X_min + d + return X + + +def Tsat_pX(p, X): + return newton(psat_TX, p, X) + + +def newton(func, p, X): + d = 1e-3 + # T = 300 + T_min = libr.trivial_keyed_output(CP.iT_min) + T_max = libr.trivial_keyed_output(CP.iT_max) + T = (T_min + T_max) / 2 + while True: + res = func(T, X) - p + deriv = (func(T + d, X) - func(T -d, X)) / (2 * d) + T -= res / deriv + + if T >= T_max: + T = T_max - d + elif T <= T_min: + T = T_min + d + + if abs(res) < 1e-6: + return T + + +# x_range = np.linspace(0, 0.75, 76) +# for x in x_range: +# psat_TX(300, x) +# Tsat_pX(1e4, x) + + +def Xsat_Tp(libr, T, p): + X_min = 0 + X_max = 0.75 + X = .3 + d = 1e-5 + while True: + res = psat_TX(libr, T, X) - p + deriv = (psat_TX(libr, T, X + d) - psat_TX(libr, T, X - d)) / (2 * d) + X -= res / deriv + # print(X) + + if abs(res) < 1e-6: + if X < X_min: + return X_min + elif X > X_max: + return X_max + break + + if X >= X_max: + X = X_max - d + elif X <= X_min: + X = X_min + d + return X + +# X = Xsat_Tp(350, 1e4) +# libr.set_mass_fractions([X]) +# print(X) + +# p_max = psat_TX(400, 0) +# print(p_max, Xsat_Tp(400, p_max)) +# for p in np.linspace(1e4, p_max): +# # print(p) +# print(p, Xsat_Tp(400, p)) + +# libr.update(CP.PT_INPUTS, 1e4, 400) + +# Tsat_pX(1e4, .3) +# libr +# T = libr.s +# libr.update(CP.QT_INPUTS, 0, T) + +# print(libr.trivial_keyed_output(CP.iT_min)) +# print(libr.trivial_keyed_output(CP.iT_freeze)) +# print(libr.trivial_keyed_output(CP.iP)) + +# watercycle +T_evap = 275.15 +p_evap = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_evap, "water") + +# T_cond = 309.15 +# p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") + +# # heat source temperature +# T_sol_des_out = 350 + +# # heat sink temperature +# T_sol_abs_out = 306.15 + +# x_rich = Xsat_Tp(T_sol_abs_out, p_evap) +# print(x_rich) +# x_water_rich = 1 - x_rich +# T_water_abs_in = T_evap +# h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") + +# h_sol_abs_out = libr.hmass() +# s_sol_abs_out = libr.smass() + +# eta_s_pump = 0.9 +# libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) +# h_pump_out_s = libr.hmass() +# h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump + +# x_poor = Xsat_Tp(T_sol_des_out, p_cond) +# h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) + +# delta_T = 0 +# T_water_des_out = T_sol_des_out - delta_T +# h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") +# print(x_poor) +# x_water_poor = 1 - x_poor + +# m_water = 1 + +# m_rich = (m_water - m_water * x_water_poor) / (x_water_rich - 1) +# m_poor = m_rich - m_water + +# delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out +# delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out + +# print(delta_H_abs) +# print(delta_H_des) + +#libr.update(CP.PT_INPUTS, p_cond, T_sol_des_out) +#print(libr.hmass()) +#print(Tsat_pX(p_cond, .5)) + +from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution + + +class Absorber(Component): + + def __init__(self, label, **kwargs): + super().__init__(label, **kwargs) + + @staticmethod + def inlets(): + return ["in1", "in2"] + + @staticmethod + def outlets(): + return ["out1"] + + def get_parameters(self): + return {} + + def get_mandatory_constraints(self): + return { + 'mass_flow_constraints': { + 'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv, + 'constant_deriv': True,# 'latex': self.mass_flow_func_doc, + 'num_eq': 1}, + 'fluid_constraints': { + 'func': self.fluid_func, 'deriv': self.fluid_deriv, + 'constant_deriv': False,# 'latex': self.fluid_func_doc, + 'num_eq': 1}, + 'pressure_constraints': { + 'func': self.pressure_equality_func, + 'deriv': self.pressure_equality_deriv, + 'constant_deriv': True, + 'latex': self.pressure_equality_func_doc, + 'num_eq': 2}, + "saturation_constraints": { + "func": self.saturated_solution_out_func, + "deriv": self.saturated_solution_out_deriv, + "constant_deriv": False, + "num_eq": 2 + } + } + + def mass_flow_func(self): + return self.inl[0].m.val_SI + self.inl[1].m.val_SI - self.outl[0].m.val_SI + + def mass_flow_deriv(self, k): + for c in self.inl: + if c.m.is_var: + self.jacobian[k, c.m.J_col] = 1 + + if self.outl[0].m.is_var: + self.jacobian[k, self.outl[0].m.J_col] = -1 + + def fluid_func(self): + return self.inl[0].m.val_SI + self.inl[1].m.val_SI * self.inl[1].fluid.val["water"] - self.outl[0].m.val_SI * self.outl[0].fluid.val["water"] + + def fluid_deriv(self, increment_filter, k): + if self.inl[0].m.is_var: + self.jacobian[k, self.inl[0].m.J_col] = 1 + if self.inl[1].m.is_var: + self.jacobian[k, self.inl[1].m.J_col] = self.inl[1].fluid.val["water"] + if self.outl[0].m.is_var: + self.jacobian[k, self.outl[0].m.J_col] = -self.outl[0].fluid.val["water"] + + + def pressure_equality_func(self): + residual = [] + for c in self.inl : + residual += [c.p.val_SI - self.outl[0].p.val_SI] + return residual + + def pressure_equality_deriv(self, k): + r""" + Calculate partial derivatives for all pressure equations. + + Returns + ------- + deriv : ndarray + Matrix with partial derivatives for the fluid equations. + """ + for c in self.inl: + if c.p.is_var: + self.jacobian[k, c.p.J_col] = 1 + if self.outl[0].p.is_var: + self.jacobian[k, self.outl[0].p.J_col] = -1 + k += 1 + + def saturated_solution_out_func(self): + outl = self.outl[0] + x_previous = outl.fluid.val["LiBr"] + T = outl.calc_T() + x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data) + return [x_libr - outl.fluid.val["LiBr"], 1 - x_libr - outl.fluid.val["water"]] + + def saturated_solution_out_deriv(self, increment_filter, k): + outl = self.outl[0] + self.jacobian[k, outl.fluid.J_col["LiBr"]] = -1 + self.jacobian[k + 1, outl.fluid.J_col["water"]] = -1 + # pass + + + @staticmethod + def is_branch_source(): + return True + + def start_branch(self): + outconn = self.outl[0] + branch = { + "connections": [outconn], + "components": [self, outconn.target], + "subbranches": {} + } + outconn.target.propagate_to_target(branch) + + return {outconn.label: branch} + + def propagate_to_target(self, branch): + return + + def propagate_wrapper_to_target(self, branch): + if branch["connections"][-1] == self.inl[0]: + return + + if self in branch["components"]: + return + + outconn = self.outl[0] + branch["connections"] += [outconn] + outconn.target.propagate_wrapper_to_target(branch) + + +from tespy.components import Source, Sink +from tespy.networks import Network +from tespy.connections import Connection + + +nw = Network() + +water = Source("water source") +rich = Sink("rich solution") +poor = Source("poor solution") + +absorber = Absorber("absorber") + +c1 = Connection(water, "out1", absorber, "in1", label="1") +c2 = Connection(poor, "out1", absorber, "in2", label="2") +c3 = Connection(absorber, "out1", rich, "in1", label="3") + +nw.add_conns(c1, c2, c3) + +c1.set_attr(fluid={"water": 1}, p=p_evap, x=1) +c2.set_attr(fluid={"INCOMP::LiBr": 0.6, "water": 0.4}, m=1, h=110000, mixing_rule="incomp-solution") +c3.set_attr(h=21000, p0=p_evap) + +nw.solve("design") +# nw.solve("design", init_only=True) + +nw + From bde0d57eb0679a308f9b86220b717d1776e958d4 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Tue, 17 Oct 2023 20:42:02 +0200 Subject: [PATCH 03/22] Add functions for saturated solution at outlet 1 --- src/tespy/tools/fluid_properties/mixtures.py | 8 +- test_libr.py | 192 +++++++++++++------ 2 files changed, 139 insertions(+), 61 deletions(-) diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index 0cdc4c240..b5cb64592 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -44,6 +44,7 @@ def xsat_pT_incomp_solution(p=None, T=None, fluid_data=None, **kwargs): res = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X) - p upper = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X + d) lower = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X - d) + fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([X]) deriv = (upper - lower) / (2 * d) X -= res / deriv @@ -122,12 +123,15 @@ def h_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): x_old = fluid_data["LiBr"]["mass_fraction"] x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs) x_water_gas = x_new - x_old - h = fluid_data["LiBr"]["wrapper"].h_pT(p, T - .005) * (1 - x_water_gas) + h = fluid_data["LiBr"]["wrapper"].h_pT(p + 1e-6, T) * (1 - x_water_gas) h += fluid_data["water"]["wrapper"].h_QT(1, T) * x_water_gas fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([x_old]) return h else: - return fluid_data["LiBr"]["wrapper"].h_pT(p, T) + try: + return fluid_data["LiBr"]["wrapper"].h_pT(p, T) + except ValueError: + return fluid_data["LiBr"]["wrapper"].h_pT(p + 1e-6, T) def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): diff --git a/test_libr.py b/test_libr.py index f10a64c20..65fa10fcd 100644 --- a/test_libr.py +++ b/test_libr.py @@ -16,20 +16,20 @@ # saturation pressure by temperature and libr mass fraction x -def psat_TX(libr, T, X): +def psat_TX(T, X): libr.set_mass_fractions([X]) libr.update(CP.QT_INPUTS, 0, T) return libr.p() -def Xsat_Tp(libr, T, p): +def Xsat_Tp(T, p): X_min = 0 X_max = 0.75 X = .3 d = 1e-5 while True: - res = psat_TX(libr, T, X) - p - deriv = (psat_TX(libr, T, X + d) - psat_TX(libr, T, X - d)) / (2 * d) + res = psat_TX(T, X) - p + deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) X -= res / deriv # print(X) @@ -77,14 +77,14 @@ def newton(func, p, X): # Tsat_pX(1e4, x) -def Xsat_Tp(libr, T, p): +def Xsat_Tp(T, p): X_min = 0 X_max = 0.75 X = .3 d = 1e-5 while True: - res = psat_TX(libr, T, X) - p - deriv = (psat_TX(libr, T, X + d) - psat_TX(libr, T, X - d)) / (2 * d) + res = psat_TX(T, X) - p + deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) X -= res / deriv # print(X) @@ -126,52 +126,46 @@ def Xsat_Tp(libr, T, p): T_evap = 275.15 p_evap = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_evap, "water") -# T_cond = 309.15 -# p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") +T_cond = 309.15 +p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") -# # heat source temperature -# T_sol_des_out = 350 +# heat source temperature +T_sol_des_out = 350 -# # heat sink temperature -# T_sol_abs_out = 306.15 +# heat sink temperature +T_sol_abs_out = 306.15 -# x_rich = Xsat_Tp(T_sol_abs_out, p_evap) -# print(x_rich) -# x_water_rich = 1 - x_rich -# T_water_abs_in = T_evap -# h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") +x_rich = Xsat_Tp(T_sol_abs_out, p_evap) -# h_sol_abs_out = libr.hmass() -# s_sol_abs_out = libr.smass() +x_water_rich = 1 - x_rich +T_water_abs_in = T_evap +h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") -# eta_s_pump = 0.9 -# libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) -# h_pump_out_s = libr.hmass() -# h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump +h_sol_abs_out = libr.hmass() +s_sol_abs_out = libr.smass() -# x_poor = Xsat_Tp(T_sol_des_out, p_cond) -# h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) +eta_s_pump = 0.9 +libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) +h_pump_out_s = libr.hmass() +h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump -# delta_T = 0 -# T_water_des_out = T_sol_des_out - delta_T -# h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") -# print(x_poor) -# x_water_poor = 1 - x_poor +x_poor = Xsat_Tp(T_sol_des_out, p_cond) +h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) -# m_water = 1 +delta_T = 0 +T_water_des_out = T_sol_des_out - delta_T +h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") -# m_rich = (m_water - m_water * x_water_poor) / (x_water_rich - 1) -# m_poor = m_rich - m_water +x_water_poor = 1 - x_poor -# delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out -# delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out +m_water = 1 +m_rich = m_water * (1 - x_water_poor) / (x_water_rich - x_water_poor) -# print(delta_H_abs) -# print(delta_H_des) +m_poor = m_rich - m_water + +delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out +delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out -#libr.update(CP.PT_INPUTS, p_cond, T_sol_des_out) -#print(libr.hmass()) -#print(Tsat_pX(p_cond, .5)) from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution @@ -208,12 +202,12 @@ def get_mandatory_constraints(self): 'constant_deriv': True, 'latex': self.pressure_equality_func_doc, 'num_eq': 2}, - "saturation_constraints": { - "func": self.saturated_solution_out_func, - "deriv": self.saturated_solution_out_deriv, - "constant_deriv": False, - "num_eq": 2 - } + # "saturation_constraints": { + # "func": self.saturated_solution_out_func, + # "deriv": self.saturated_solution_out_deriv, + # "constant_deriv": False, + # "num_eq": 2 + # } } def mass_flow_func(self): @@ -261,19 +255,36 @@ def pressure_equality_deriv(self, k): self.jacobian[k, self.outl[0].p.J_col] = -1 k += 1 - def saturated_solution_out_func(self): + def saturated_solution_water_func(self): + outl = self.outl[0] + return 1 - outl.fluid.val["LiBr"] - outl.fluid.val["water"] + + def saturated_solution_water_deriv(self, increment_filter, k): + outl = self.outl[0] + if "water" in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col["water"]] = -1 + if "LiBr" in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col["LiBr"]] = -1 + # pass + + def saturated_solution_libr_func(self): outl = self.outl[0] x_previous = outl.fluid.val["LiBr"] T = outl.calc_T() x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data) - return [x_libr - outl.fluid.val["LiBr"], 1 - x_libr - outl.fluid.val["water"]] + outl.fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([x_previous]) + return x_libr - outl.fluid.val["LiBr"] - def saturated_solution_out_deriv(self, increment_filter, k): + def saturated_solution_libr_deriv(self, increment_filter, k): outl = self.outl[0] - self.jacobian[k, outl.fluid.J_col["LiBr"]] = -1 - self.jacobian[k + 1, outl.fluid.J_col["water"]] = -1 - # pass - + if outl.p.is_var: + deriv = self.numeric_deriv(self.saturated_solution_libr_func, "p", outl) + self.jacobian[k, outl.p.J_col] = deriv + if outl.h.is_var: + deriv = self.numeric_deriv(self.saturated_solution_libr_func, "h", outl) + self.jacobian[k, outl.h.J_col] = deriv + if "LiBr" in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col["LiBr"]] = -1 @staticmethod def is_branch_source(): @@ -304,6 +315,38 @@ def propagate_wrapper_to_target(self, branch): branch["connections"] += [outconn] outconn.target.propagate_wrapper_to_target(branch) + def new_constraints(self): + constraints = { + 'mass_flow_constraints': { + 'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv, + 'constant_deriv': True,# 'latex': self.mass_flow_func_doc, + 'num_eq': 1}, + 'fluid_constraints': { + 'func': self.fluid_func, 'deriv': self.fluid_deriv, + 'constant_deriv': False,# 'latex': self.fluid_func_doc, + 'num_eq': 1}, + 'pressure_constraints': { + 'func': self.pressure_equality_func, + 'deriv': self.pressure_equality_deriv, + 'constant_deriv': True, + 'latex': self.pressure_equality_func_doc, + 'num_eq': 2}, + "saturation_constraints_libr": { + "func": self.saturated_solution_libr_func, + "deriv": self.saturated_solution_libr_deriv, + "constant_deriv": False, + "num_eq": 1 + }, + } + if "LiBr" in self.outl[0].fluid.is_var: + constraints["saturation_constraints_water"] = { + "func": self.saturated_solution_water_func, + "deriv": self.saturated_solution_water_deriv, + "constant_deriv": False, + "num_eq": 1 + } + return constraints + from tespy.components import Source, Sink from tespy.networks import Network @@ -324,12 +367,43 @@ def propagate_wrapper_to_target(self, branch): nw.add_conns(c1, c2, c3) -c1.set_attr(fluid={"water": 1}, p=p_evap, x=1) -c2.set_attr(fluid={"INCOMP::LiBr": 0.6, "water": 0.4}, m=1, h=110000, mixing_rule="incomp-solution") -c3.set_attr(h=21000, p0=p_evap) +c1.set_attr(fluid={"water": 1}, p=p_evap, m=1, x=1) +c2.set_attr(fluid={"water": x_water_poor}, h=h_poor, mixing_rule="incomp-solution") +c3.set_attr(fluid={"INCOMP::LiBr": x_rich}, h=h_sol_abs_out, p0=p_evap) +nw.set_attr() +nw.solve("design") +# how does the reference temperature affect the energy balance +p_ref = 1e5 +for T_ref in range(274, 400): + h_water_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, "water") + h_poor_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_poor}]") + h_rich_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_rich}]") + print((c1.h.val_SI - h_water_ref) * c1.m.val_SI + (c2.h.val_SI - h_poor_ref) * c2.m.val_SI - (c3.h.val_SI - h_rich_ref) * c3.m.val_SI) + +# some checks +assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) +assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) + +# replace initial draft of constraints -> maybe add saturation constraint as parameter +absorber.get_mandatory_constraints = absorber.new_constraints +# now there is two additional equations in case both fluids at outlet 1 are variable +# and one additional equation if none of the fluids is variale +# fluid composition as function of saturation temperature and pressure +c3.fluid.is_set = set() +c1.set_attr(p=None) +c3.set_attr(h=None, T=306.15, p=p_evap) +c3.h.val0 = c3.h.val_SI nw.solve("design") -# nw.solve("design", init_only=True) +# check results, should be identical to previous ones +assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) +assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) -nw +# fix the either the water or the libr mass fraction and the temperature or pressure +# gives us the other respective value, pressure in this case +c3.set_attr(fluid={"INCOMP::LiBr": 0.5}, T=295, p=None) +nw.solve("design") +assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) +assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) +print(f"T3: {c3.T.val_SI}, p3: {c3.p.val_SI}, fluid3: {c3.fluid.val}") From 7c713793cfa7226553eaac1d1f5e9ba610e55fe9 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Tue, 17 Oct 2023 20:42:45 +0200 Subject: [PATCH 04/22] Replace hard coded fluid names --- test_libr.py | 46 ++++++++++++++++++++++++---------------------- 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/test_libr.py b/test_libr.py index 65fa10fcd..d7d8c0345 100644 --- a/test_libr.py +++ b/test_libr.py @@ -183,6 +183,11 @@ def inlets(): def outlets(): return ["out1"] + def preprocess(self, num_vars): + self.h2o = "water" + self.solvent = "LiBr" + super().preprocess(num_vars) + def get_parameters(self): return {} @@ -202,12 +207,6 @@ def get_mandatory_constraints(self): 'constant_deriv': True, 'latex': self.pressure_equality_func_doc, 'num_eq': 2}, - # "saturation_constraints": { - # "func": self.saturated_solution_out_func, - # "deriv": self.saturated_solution_out_deriv, - # "constant_deriv": False, - # "num_eq": 2 - # } } def mass_flow_func(self): @@ -222,16 +221,19 @@ def mass_flow_deriv(self, k): self.jacobian[k, self.outl[0].m.J_col] = -1 def fluid_func(self): - return self.inl[0].m.val_SI + self.inl[1].m.val_SI * self.inl[1].fluid.val["water"] - self.outl[0].m.val_SI * self.outl[0].fluid.val["water"] + return ( + self.inl[0].m.val_SI + + self.inl[1].m.val_SI * self.inl[1].fluid.val[self.h2o] + - self.outl[0].m.val_SI * self.outl[0].fluid.val[self.h2o] + ) def fluid_deriv(self, increment_filter, k): if self.inl[0].m.is_var: self.jacobian[k, self.inl[0].m.J_col] = 1 if self.inl[1].m.is_var: - self.jacobian[k, self.inl[1].m.J_col] = self.inl[1].fluid.val["water"] + self.jacobian[k, self.inl[1].m.J_col] = self.inl[1].fluid.val[self.h2o] if self.outl[0].m.is_var: - self.jacobian[k, self.outl[0].m.J_col] = -self.outl[0].fluid.val["water"] - + self.jacobian[k, self.outl[0].m.J_col] = -self.outl[0].fluid.val[self.h2o] def pressure_equality_func(self): residual = [] @@ -257,23 +259,23 @@ def pressure_equality_deriv(self, k): def saturated_solution_water_func(self): outl = self.outl[0] - return 1 - outl.fluid.val["LiBr"] - outl.fluid.val["water"] + return 1 - outl.fluid.val[self.solvent] - outl.fluid.val[self.h2o] def saturated_solution_water_deriv(self, increment_filter, k): outl = self.outl[0] - if "water" in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col["water"]] = -1 - if "LiBr" in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col["LiBr"]] = -1 + if self.h2o in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[self.h2o]] = -1 + if self.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[self.solvent]] = -1 # pass def saturated_solution_libr_func(self): outl = self.outl[0] - x_previous = outl.fluid.val["LiBr"] - T = outl.calc_T() - x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data) - outl.fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([x_previous]) - return x_libr - outl.fluid.val["LiBr"] + x_previous = outl.fluid.val[self.solvent] + T = outl.calc_T(solvent=self.solvent) + x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=self.solvent) + outl.fluid_data[self.solvent]["wrapper"].AS.set_mass_fractions([x_previous]) + return x_libr - outl.fluid.val[self.solvent] def saturated_solution_libr_deriv(self, increment_filter, k): outl = self.outl[0] @@ -283,8 +285,8 @@ def saturated_solution_libr_deriv(self, increment_filter, k): if outl.h.is_var: deriv = self.numeric_deriv(self.saturated_solution_libr_func, "h", outl) self.jacobian[k, outl.h.J_col] = deriv - if "LiBr" in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col["LiBr"]] = -1 + if self.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[self.solvent]] = -1 @staticmethod def is_branch_source(): From 4b8cc117db18b8a6627282b8fd99b02e550fcfc9 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Fri, 20 Oct 2023 15:11:24 +0200 Subject: [PATCH 05/22] Add a solvent to incompressible solution mixtures --- src/tespy/connections/connection.py | 21 ++++++++-------- src/tespy/networks/network.py | 19 +++++++++++++-- src/tespy/tools/fluid_properties/functions.py | 24 +++++++++---------- src/tespy/tools/fluid_properties/helpers.py | 5 ++-- src/tespy/tools/fluid_properties/mixtures.py | 22 +++++++++-------- 5 files changed, 55 insertions(+), 36 deletions(-) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index 5089f06bf..ec47dc74a 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -271,6 +271,7 @@ def __init__(self, source, outlet_id, target, inlet_id, self.property_data0 = [x + '0' for x in self.property_data.keys()] self.__dict__.update(self.property_data) self.mixing_rule = None + self.solvent = None msg = ( f"Created connection from {self.source.label} ({self.source_id}) " f"to {self.target.label} ({self.target_id})." @@ -457,8 +458,8 @@ def set_attr(self, **kwargs): else: self.__dict__.update({key: kwargs[key]}) - elif key == "mixing_rule": - self.mixing_rule = kwargs[key] + elif key in ["mixing_rule", "solvent"]: + self.__dict__.update({key: kwargs[key]}) # invalid keyword else: @@ -720,7 +721,7 @@ def build_fluid_data(self): } for fluid in self.fluid.val } if self.mixing_rule == "incomp-solution": - self.fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([self.fluid.val["LiBr"]]) + self.fluid_data[self.solvent]["wrapper"].AS.set_mass_fractions([self.fluid.val[self.solvent]]) def primary_ref_func(self, k, **kwargs): variable = kwargs["variable"] @@ -741,7 +742,7 @@ def primary_ref_deriv(self, k, **kwargs): self.jacobian[k, ref.obj.get_attr(variable).J_col] = -ref.factor def calc_T(self, T0=None): - return T_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0) + return T_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0, solvent=self.solvent) def T_func(self, k, **kwargs): self.residual[k] = self.calc_T() - self.T.val_SI @@ -749,15 +750,15 @@ def T_func(self, k, **kwargs): def T_deriv(self, k, **kwargs): if self.p.is_var: self.jacobian[k, self.p.J_col] = ( - dT_mix_dph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI) + dT_mix_dph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI, solvent=self.solvent) ) if self.h.is_var: self.jacobian[k, self.h.J_col] = ( - dT_mix_pdh(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI) + dT_mix_pdh(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, self.T.val_SI, solvent=self.solvent) ) for fluid in self.fluid.is_var: self.jacobian[k, self.fluid.J_col[fluid]] = dT_mix_ph_dfluid( - self.p.val_SI, self.h.val_SI, fluid, self.fluid_data, self.mixing_rule + self.p.val_SI, self.h.val_SI, fluid, self.fluid_data, self.mixing_rule, solvent=self.solvent ) def T_ref_func(self, k, **kwargs): @@ -772,16 +773,16 @@ def T_ref_deriv(self, k, **kwargs): ref = self.T_ref.ref if ref.obj.p.is_var: self.jacobian[k, ref.obj.p.J_col] = -( - dT_mix_dph(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule) + dT_mix_dph(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule, solvent=ref.obj.solvent) ) * ref.factor if ref.obj.h.is_var: self.jacobian[k, ref.obj.h.J_col] = -( - dT_mix_pdh(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule) + dT_mix_pdh(ref.obj.p.val_SI, ref.obj.h.val_SI, ref.obj.fluid_data, ref.obj.mixing_rule, solvent=ref.obj.solvent) ) * ref.factor for fluid in ref.obj.fluid.is_var: if not self._increment_filter[ref.obj.fluid.J_col[fluid]]: self.jacobian[k, ref.obj.fluid.J_col[fluid]] = -dT_mix_ph_dfluid( - ref.obj.p.val_SI, ref.obj.h.val_SI, fluid, ref.obj.fluid_data, ref.obj.mixing_rule + ref.obj.p.val_SI, ref.obj.h.val_SI, fluid, ref.obj.fluid_data, ref.obj.mixing_rule, solvent=ref.obj.solvent ) def calc_viscosity(self, T0=None): diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index 405c0b2eb..bf6f66521 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -893,16 +893,30 @@ def propagate_fluid_wrappers(self): if f in c.fluid.back_end: back_ends[f] = c.fluid.back_end[f] + solvents = [ + c.solvent for c in all_connections + if c.solvent is not None + ] mixing_rules = [ c.mixing_rule for c in all_connections if c.mixing_rule is not None ] mixing_rule = set(mixing_rules) + solvent = set(solvents) + if len(solvent) > 1: + msg = "You have provided more than one solvent." + raise hlp.TESPyNetworkError(msg) + elif len(solvent) == 0: + solvent = None + else: + solvent = solvent.pop() if len(mixing_rule) > 1: msg = "You have provided more than one mixing rule." raise hlp.TESPyNetworkError(msg) elif len(mixing_rule) == 0: - mixing_rule = set(["ideal-cond"]) + mixing_rule = "ideal-cond" + else: + mixing_rule = mixing_rule.pop() if not any_fluids_set: msg = "You are missing fluid specifications." @@ -920,7 +934,8 @@ def propagate_fluid_wrappers(self): raise hlp.TESPyNetworkError(msg) for c in all_connections: - c.mixing_rule = list(mixing_rule)[0] + c.solvent = solvent + c.mixing_rule = mixing_rule c._potential_fluids = potential_fluids if num_potential_fluids == 1: f = list(potential_fluids)[0] diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py index f97ce42d1..1845107ed 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/src/tespy/tools/fluid_properties/functions.py @@ -85,39 +85,39 @@ def calc_chemical_exergy(pamb, Tamb, fluid_data, Chem_Ex, mixing_rule=None, T0=N return EXERGY_CHEMICAL[mixing_rule](pamb, Tamb, fluid_data, Chem_Ex) -def T_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): +def T_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].T_ph(p, h) else: _check_mixing_rule(mixing_rule, T_MIX_PH_REVERSE, "temperature (from enthalpy)") - kwargs = { + kwargs.update({ "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): +def dT_mix_pdh(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs): d = 1e-1 - upper = T_mix_ph(p, h + d, fluid_data, mixing_rule=mixing_rule, T0=T0) - lower = T_mix_ph(p, h - d, fluid_data, mixing_rule=mixing_rule, T0=upper) + upper = T_mix_ph(p, h + d, fluid_data, mixing_rule=mixing_rule, T0=T0, **kwargs) + lower = T_mix_ph(p, h - d, fluid_data, mixing_rule=mixing_rule, T0=upper, **kwargs) return (upper - lower) / (2 * d) -def dT_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None): +def dT_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs): d = 1e-1 - upper = T_mix_ph(p + d, h, fluid_data, mixing_rule=mixing_rule, T0=T0) - lower = T_mix_ph(p - d, h, fluid_data, mixing_rule=mixing_rule, T0=upper) + upper = T_mix_ph(p + d, h, fluid_data, mixing_rule=mixing_rule, T0=T0, **kwargs) + lower = T_mix_ph(p - d, h, fluid_data, mixing_rule=mixing_rule, T0=upper, **kwargs) return (upper - lower) / (2 * d) -def dT_mix_ph_dfluid(p, h, fluid, fluid_data, mixing_rule=None, T0=None): +def dT_mix_ph_dfluid(p, h, fluid, fluid_data, mixing_rule=None, T0=None, **kwargs): d = 1e-5 fluid_data[fluid]["mass_fraction"] += d - upper = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=T0) + upper = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=T0, **kwargs) fluid_data[fluid]["mass_fraction"] -= 2 * d - lower = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=upper) + lower = T_mix_ph(p, h, fluid_data, mixing_rule=mixing_rule, T0=upper, **kwargs) fluid_data[fluid]["mass_fraction"] += d return (upper - lower) / (2 * d) diff --git a/src/tespy/tools/fluid_properties/helpers.py b/src/tespy/tools/fluid_properties/helpers.py index 226648579..3b6f34e55 100644 --- a/src/tespy/tools/fluid_properties/helpers.py +++ b/src/tespy/tools/fluid_properties/helpers.py @@ -64,7 +64,7 @@ def get_molar_fractions(fluid_data): return {key: value / molarflow_sum for key, value in molarflow.items()} -def inverse_temperature_mixture(p=None, target_value=None, fluid_data=None, T0=None, f=None): +def inverse_temperature_mixture(p=None, target_value=None, fluid_data=None, T0=None, f=None, **kwargs): # calculate the fluid properties for fluid mixtures valmin, valmax = get_mixture_temperature_range(fluid_data) if T0 is None: @@ -74,13 +74,14 @@ def inverse_temperature_mixture(p=None, target_value=None, fluid_data=None, T0=N "p": p, "fluid_data": fluid_data, "T": T0, "function": f, "parameter": "T" , "delta": 0.01 } + function_kwargs.update(**kwargs) return newton_with_kwargs( central_difference, target_value, val0=T0, valmin=valmin, valmax=valmax, - **function_kwargs + **function_kwargs, ) diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index b5cb64592..7c3d15d41 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -40,11 +40,12 @@ def xsat_pT_incomp_solution(p=None, T=None, fluid_data=None, **kwargs): X = .3 d = 1e-5 iter = 0 + solvent = kwargs["solvent"] while True: - res = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X) - p - upper = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X + d) - lower = fluid_data["LiBr"]["wrapper"].psat_Tx(T, X - d) - fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([X]) + res = fluid_data[solvent]["wrapper"].psat_Tx(T, X) - p + upper = fluid_data[solvent]["wrapper"].psat_Tx(T, X + d) + lower = fluid_data[solvent]["wrapper"].psat_Tx(T, X - d) + fluid_data[solvent]["wrapper"].AS.set_mass_fractions([X]) deriv = (upper - lower) / (2 * d) X -= res / deriv @@ -119,19 +120,20 @@ def h_mix_pT_incompressible(p, T, fluid_data, **kwargs): def h_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): - if p < fluid_data["LiBr"]["wrapper"].p_sat(T): - x_old = fluid_data["LiBr"]["mass_fraction"] + solvent = kwargs["solvent"] + if p < fluid_data[solvent]["wrapper"].p_sat(T): + x_old = fluid_data[solvent]["mass_fraction"] x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs) x_water_gas = x_new - x_old - h = fluid_data["LiBr"]["wrapper"].h_pT(p + 1e-6, T) * (1 - x_water_gas) + h = fluid_data[solvent]["wrapper"].h_pT(p + 1e-6, T) * (1 - x_water_gas) h += fluid_data["water"]["wrapper"].h_QT(1, T) * x_water_gas - fluid_data["LiBr"]["wrapper"].AS.set_mass_fractions([x_old]) + fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old]) return h else: try: - return fluid_data["LiBr"]["wrapper"].h_pT(p, T) + return fluid_data[solvent]["wrapper"].h_pT(p, T) except ValueError: - return fluid_data["LiBr"]["wrapper"].h_pT(p + 1e-6, T) + return fluid_data[solvent]["wrapper"].h_pT(p + 1e-6, T) def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): From 01319ce2b2a05b3c4f19f0b0d28ce01b58912a54 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Fri, 20 Oct 2023 15:13:48 +0200 Subject: [PATCH 06/22] Call the solvent from the connection --- test_libr.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/test_libr.py b/test_libr.py index d7d8c0345..1774306e1 100644 --- a/test_libr.py +++ b/test_libr.py @@ -185,7 +185,6 @@ def outlets(): def preprocess(self, num_vars): self.h2o = "water" - self.solvent = "LiBr" super().preprocess(num_vars) def get_parameters(self): @@ -259,23 +258,24 @@ def pressure_equality_deriv(self, k): def saturated_solution_water_func(self): outl = self.outl[0] - return 1 - outl.fluid.val[self.solvent] - outl.fluid.val[self.h2o] + return 1 - outl.fluid.val[outl.solvent] - outl.fluid.val[self.h2o] + def saturated_solution_water_deriv(self, increment_filter, k): outl = self.outl[0] if self.h2o in outl.fluid.is_var: self.jacobian[k, outl.fluid.J_col[self.h2o]] = -1 - if self.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[self.solvent]] = -1 + if outl.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -1 # pass def saturated_solution_libr_func(self): outl = self.outl[0] - x_previous = outl.fluid.val[self.solvent] - T = outl.calc_T(solvent=self.solvent) - x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=self.solvent) - outl.fluid_data[self.solvent]["wrapper"].AS.set_mass_fractions([x_previous]) - return x_libr - outl.fluid.val[self.solvent] + x_previous = outl.fluid.val[outl.solvent] + T = outl.calc_T() + x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=outl.solvent) + outl.fluid_data[outl.solvent]["wrapper"].AS.set_mass_fractions([x_previous]) + return x_libr - outl.fluid.val[outl.solvent] def saturated_solution_libr_deriv(self, increment_filter, k): outl = self.outl[0] @@ -285,8 +285,8 @@ def saturated_solution_libr_deriv(self, increment_filter, k): if outl.h.is_var: deriv = self.numeric_deriv(self.saturated_solution_libr_func, "h", outl) self.jacobian[k, outl.h.J_col] = deriv - if self.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[self.solvent]] = -1 + if outl.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -1 @staticmethod def is_branch_source(): @@ -370,7 +370,7 @@ def new_constraints(self): nw.add_conns(c1, c2, c3) c1.set_attr(fluid={"water": 1}, p=p_evap, m=1, x=1) -c2.set_attr(fluid={"water": x_water_poor}, h=h_poor, mixing_rule="incomp-solution") +c2.set_attr(fluid={"water": x_water_poor}, h=h_poor, mixing_rule="incomp-solution", solvent="LiBr") c3.set_attr(fluid={"INCOMP::LiBr": x_rich}, h=h_sol_abs_out, p0=p_evap) nw.set_attr() nw.solve("design") @@ -403,7 +403,7 @@ def new_constraints(self): # fix the either the water or the libr mass fraction and the temperature or pressure # gives us the other respective value, pressure in this case -c3.set_attr(fluid={"INCOMP::LiBr": 0.5}, T=295, p=None) +c3.set_attr(fluid={"INCOMP::LiBr": 0.45}, T=295, p=None) nw.solve("design") assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) From f77b1c5e4ceade5896d2b8325751eea16613774c Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Fri, 20 Oct 2023 19:02:51 +0200 Subject: [PATCH 07/22] Add **kwargs for fluid property functions to make solvent usable --- src/tespy/components/turbomachinery/pump.py | 5 +-- src/tespy/connections/connection.py | 4 +-- src/tespy/tools/fluid_properties/functions.py | 32 +++++++++---------- src/tespy/tools/fluid_properties/mixtures.py | 24 ++++++++++++-- test_libr.py | 16 +++++++++- 5 files changed, 58 insertions(+), 23 deletions(-) diff --git a/src/tespy/components/turbomachinery/pump.py b/src/tespy/components/turbomachinery/pump.py index 2fd0568f8..32d68766e 100644 --- a/src/tespy/components/turbomachinery/pump.py +++ b/src/tespy/components/turbomachinery/pump.py @@ -211,7 +211,8 @@ def eta_s_func(self): o.p.val_SI, i.fluid_data, i.mixing_rule, - T0=None + T0=None, + solvent=i.solvent ) - self.inl[0].h.val_SI ) ) @@ -511,7 +512,7 @@ def calc_parameters(self): o.p.val_SI, i.fluid_data, i.mixing_rule, - T0=None + T0=None, solvent=i.solvent ) - self.inl[0].h.val_SI ) / (o.h.val_SI - i.h.val_SI) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index ec47dc74a..bb495716c 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -787,14 +787,14 @@ def T_ref_deriv(self, k, **kwargs): def calc_viscosity(self, T0=None): try: - return viscosity_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0) + return viscosity_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0, solvent=self.solvent) except NotImplementedError: return np.nan def calc_vol(self, T0=None): try: - return v_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0) + return v_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=T0, solvent=self.solvent) except NotImplementedError: return np.nan diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py index 1845107ed..8abb4bf35 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/src/tespy/tools/fluid_properties/functions.py @@ -24,14 +24,14 @@ from .mixtures import VISCOSITY_MIX_PT_DIRECT -def isentropic(p_1, h_1, p_2, fluid_data, mixing_rule=None, T0=None): +def isentropic(p_1, h_1, p_2, fluid_data, mixing_rule=None, T0=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].isentropic(p_1, h_1, p_2) else: - s_1 = s_mix_ph(p_1, h_1, fluid_data, mixing_rule) - T_2 = T_mix_ps(p_2, s_1, fluid_data, mixing_rule) - return h_mix_pT(p_2, T_2, fluid_data, mixing_rule) + s_1 = s_mix_ph(p_1, h_1, fluid_data, mixing_rule, **kwargs) + T_2 = T_mix_ps(p_2, s_1, fluid_data, mixing_rule, **kwargs) + return h_mix_pT(p_2, T_2, fluid_data, mixing_rule, **kwargs) def calc_physical_exergy(h, s, p, pamb, Tamb, fluid_data, mixing_rule=None, T0=None): @@ -122,13 +122,13 @@ def dT_mix_ph_dfluid(p, h, fluid, fluid_data, mixing_rule=None, T0=None, **kwarg return (upper - lower) / (2 * d) -def h_mix_pT(p, T, fluid_data, mixing_rule=None): +def h_mix_pT(p, T, fluid_data, mixing_rule=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].h_pT(p, T) else: _check_mixing_rule(mixing_rule, H_MIX_PT_DIRECT, "enthalpy") - return H_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data) + return H_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data, **kwargs) def h_mix_pQ(p, Q, fluid_data, mixing_rule=None): @@ -181,39 +181,39 @@ def dT_sat_dp(p, fluid_data, mixing_rule=None): return (upper - lower) / (2 * d) -def s_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): +def s_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].s_ph(p, h) else: - T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) - return s_mix_pT(p, T, fluid_data, mixing_rule) + T = T_mix_ph(p, h , fluid_data, mixing_rule, T0, **kwargs) + return s_mix_pT(p, T, fluid_data, mixing_rule, **kwargs) -def s_mix_pT(p, T, fluid_data, mixing_rule=None): +def s_mix_pT(p, T, fluid_data, mixing_rule=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].s_pT(p, T) else: _check_mixing_rule(mixing_rule, S_MIX_PT_DIRECT, "entropy") - return S_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data) + return S_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data, **kwargs) -def T_mix_ps(p, s, fluid_data, mixing_rule=None, T0=None): +def T_mix_ps(p, s, fluid_data, mixing_rule=None, T0=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: 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 = { + kwargs.update({ "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): +def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: pure_fluid = get_pure_fluid(fluid_data) return 1 / pure_fluid["wrapper"].d_ph(p, h) @@ -236,7 +236,7 @@ def dv_mix_pdh(p, h, fluid_data, mixing_rule=None, T0=None): return (upper - lower) / (2 * d) -def v_mix_pT(p, T, fluid_data, mixing_rule=None): +def v_mix_pT(p, T, fluid_data, mixing_rule=None, **kwargs): if get_number_of_fluids(fluid_data) == 1: pure_fluid = get_pure_fluid(fluid_data) return 1 / pure_fluid["wrapper"].d_pT(p, T) diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index 7c3d15d41..c201f94f7 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -182,6 +182,23 @@ def s_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return s +def s_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): + solvent = kwargs["solvent"] + if p < fluid_data[solvent]["wrapper"].p_sat(T): + x_old = fluid_data[solvent]["mass_fraction"] + x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs) + x_water_gas = x_new - x_old + s = fluid_data[solvent]["wrapper"].s_pT(p + 1e-6, T) * (1 - x_water_gas) + s += fluid_data["water"]["wrapper"].s_QT(1, T) * x_water_gas + fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old]) + return s + else: + try: + return fluid_data[solvent]["wrapper"].s_pT(p, T) + except ValueError: + return fluid_data[solvent]["wrapper"].s_pT(p + 1e-6, T) + + def v_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): molar_fractions = get_molar_fractions(fluid_data) @@ -395,7 +412,8 @@ def cond_check(p, T, fluid_data, water_alias): 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, + "incomp-solution": s_mix_pT_incompressible_solution } @@ -403,6 +421,7 @@ 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, + "incomp-solution": h_mix_pT_incompressible_solution, "forced-gas": h_mix_pT_forced_gas } @@ -410,7 +429,8 @@ def cond_check(p, T, fluid_data, water_alias): 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, + "incomp-solution": s_mix_pT_incompressible_solution } diff --git a/test_libr.py b/test_libr.py index 1774306e1..32fd21c6b 100644 --- a/test_libr.py +++ b/test_libr.py @@ -191,7 +191,7 @@ def get_parameters(self): return {} def get_mandatory_constraints(self): - return { + constraints = { 'mass_flow_constraints': { 'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv, 'constant_deriv': True,# 'latex': self.mass_flow_func_doc, @@ -206,7 +206,21 @@ def get_mandatory_constraints(self): 'constant_deriv': True, 'latex': self.pressure_equality_func_doc, 'num_eq': 2}, + "saturation_constraints_libr": { + "func": self.saturated_solution_libr_func, + "deriv": self.saturated_solution_libr_deriv, + "constant_deriv": False, + "num_eq": 1 + }, } + if "LiBr" in self.outl[0].fluid.is_var: + constraints["saturation_constraints_water"] = { + "func": self.saturated_solution_water_func, + "deriv": self.saturated_solution_water_deriv, + "constant_deriv": False, + "num_eq": 1 + } + return constraints def mass_flow_func(self): return self.inl[0].m.val_SI + self.inl[1].m.val_SI - self.outl[0].m.val_SI From 3cfcacba813bda73899dc608bd68b488b0f0b251 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Fri, 20 Oct 2023 19:03:13 +0200 Subject: [PATCH 08/22] Calculate mixture temperatures, relax on convergence conditions for now --- src/tespy/connections/connection.py | 78 ++++++++++++++--------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index bb495716c..eb9d8db19 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -906,52 +906,52 @@ def solve(self, increment_filter): def calc_results(self): _converged = True - try: - self.T.val_SI = self.calc_T() - number_fluids = get_number_of_fluids(self.fluid_data) - if number_fluids > 1: - h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule) - if abs(h_from_T - self.h.val_SI) > ERR ** .5: - self.T.val_SI = np.nan - self.vol.val_SI = np.nan - self.v.val_SI = np.nan - self.s.val_SI = np.nan - msg = ( - "Could not find a feasible value for mixture temperature at " - f"connection {self.label}. The values for temperature, " - "specific volume, volumetric flow and entropy are set to nan." - ) - logger.error(msg) - _converged = False + self.T.val_SI = self.calc_T() + number_fluids = get_number_of_fluids(self.fluid_data) + if number_fluids > 1: + h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule, solvent=self.solvent) + print(h_from_T, self.h.val_SI) + if abs(h_from_T - self.h.val_SI) > 1e-1: + self.T.val_SI = np.nan + self.vol.val_SI = np.nan + self.v.val_SI = np.nan + self.s.val_SI = np.nan + msg = ( + "Could not find a feasible value for mixture temperature at " + f"connection {self.label}. The values for temperature, " + "specific volume, volumetric flow and entropy are set to nan." + ) + logger.error(msg) + _converged = False - else: - try: - if not self.x.is_set: - self.x.val_SI = self.calc_x() - except ValueError: - self.x.val_SI = np.nan - try: - if not self.Td_bp.is_set: - self.Td_bp.val_SI = self.calc_Td_bp() - except ValueError: - self.x.val_SI = np.nan + else: + try: + if not self.x.is_set: + self.x.val_SI = self.calc_x() + except ValueError: + self.x.val_SI = np.nan + try: + if not self.Td_bp.is_set: + self.Td_bp.val_SI = self.calc_Td_bp() + except ValueError: + self.x.val_SI = np.nan + try: if _converged: self.vol.val_SI = self.calc_vol() self.v.val_SI = self.vol.val_SI * self.m.val_SI self.s.val_SI = self.calc_s() - - for prop in fpd.keys(): - self.get_attr(prop).val = convert_from_SI( - prop, self.get_attr(prop).val_SI, self.get_attr(prop).unit - ) - - self.m.val0 = self.m.val - self.p.val0 = self.p.val - self.h.val0 = self.h.val - self.fluid.val0 = self.fluid.val.copy() - except (ValueError, KeyError): + except KeyError: pass + for prop in fpd.keys(): + self.get_attr(prop).val = convert_from_SI( + prop, self.get_attr(prop).val_SI, self.get_attr(prop).unit + ) + + self.m.val0 = self.m.val + self.p.val0 = self.p.val + self.h.val0 = self.h.val + self.fluid.val0 = self.fluid.val.copy() def check_pressure_bounds(self, fluid): if self.p.val_SI > self.fluid.wrapper[fluid]._p_max: From f94294ae3f01006324dcf5a17b3cc0e487909902 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Fri, 20 Oct 2023 19:03:28 +0200 Subject: [PATCH 09/22] Add a pump, heat exchanger and a valve to the mix --- test_libr.py | 110 +++++++++++++++++++++++++++------------------------ 1 file changed, 59 insertions(+), 51 deletions(-) diff --git a/test_libr.py b/test_libr.py index 32fd21c6b..5392fe817 100644 --- a/test_libr.py +++ b/test_libr.py @@ -331,40 +331,8 @@ def propagate_wrapper_to_target(self, branch): branch["connections"] += [outconn] outconn.target.propagate_wrapper_to_target(branch) - def new_constraints(self): - constraints = { - 'mass_flow_constraints': { - 'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv, - 'constant_deriv': True,# 'latex': self.mass_flow_func_doc, - 'num_eq': 1}, - 'fluid_constraints': { - 'func': self.fluid_func, 'deriv': self.fluid_deriv, - 'constant_deriv': False,# 'latex': self.fluid_func_doc, - 'num_eq': 1}, - 'pressure_constraints': { - 'func': self.pressure_equality_func, - 'deriv': self.pressure_equality_deriv, - 'constant_deriv': True, - 'latex': self.pressure_equality_func_doc, - 'num_eq': 2}, - "saturation_constraints_libr": { - "func": self.saturated_solution_libr_func, - "deriv": self.saturated_solution_libr_deriv, - "constant_deriv": False, - "num_eq": 1 - }, - } - if "LiBr" in self.outl[0].fluid.is_var: - constraints["saturation_constraints_water"] = { - "func": self.saturated_solution_water_func, - "deriv": self.saturated_solution_water_deriv, - "constant_deriv": False, - "num_eq": 1 - } - return constraints - -from tespy.components import Source, Sink +from tespy.components import Source, Sink, Pump, Valve, HeatExchanger from tespy.networks import Network from tespy.connections import Connection @@ -375,34 +343,45 @@ def new_constraints(self): rich = Sink("rich solution") poor = Source("poor solution") +pu = Pump("pump") +va = Valve("valve") + absorber = Absorber("absorber") c1 = Connection(water, "out1", absorber, "in1", label="1") -c2 = Connection(poor, "out1", absorber, "in2", label="2") -c3 = Connection(absorber, "out1", rich, "in1", label="3") +c2 = Connection(va, "out1", absorber, "in2", label="2") +c3 = Connection(absorber, "out1", pu, "in1", label="3") +c4 = Connection(pu, "out1", rich, "in1", label="4") +c7 = Connection(poor, "out1", va, "in1", label="7") -nw.add_conns(c1, c2, c3) +nw.add_conns(c1, c2, c3, c4, c7) + +# pu.set_attr(eta_s=0.8) c1.set_attr(fluid={"water": 1}, p=p_evap, m=1, x=1) -c2.set_attr(fluid={"water": x_water_poor}, h=h_poor, mixing_rule="incomp-solution", solvent="LiBr") -c3.set_attr(fluid={"INCOMP::LiBr": x_rich}, h=h_sol_abs_out, p0=p_evap) +c2.set_attr(fluid={"water": x_water_poor, "INCOMP::LiBr": 1 - x_water_poor}, h=h_poor, mixing_rule="incomp-solution", solvent="LiBr") +c3.set_attr(h=h_sol_abs_out, p0=p_evap) +print(h_sol_abs_out) +c4.set_attr(p=p_cond, h=37410) +c7.set_attr(p=p_cond) nw.set_attr() nw.solve("design") # how does the reference temperature affect the energy balance -p_ref = 1e5 -for T_ref in range(274, 400): - h_water_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, "water") - h_poor_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_poor}]") - h_rich_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_rich}]") - print((c1.h.val_SI - h_water_ref) * c1.m.val_SI + (c2.h.val_SI - h_poor_ref) * c2.m.val_SI - (c3.h.val_SI - h_rich_ref) * c3.m.val_SI) +# p_ref = 1e5 +# for T_ref in range(274, 400): +# h_water_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, "water") +# h_poor_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_poor}]") +# h_rich_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_rich}]") +# print((c1.h.val_SI - h_water_ref) * c1.m.val_SI + (c2.h.val_SI - h_poor_ref) * c2.m.val_SI - (c3.h.val_SI - h_rich_ref) * c3.m.val_SI) # some checks -assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) -assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) +nw.print_results() +assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 2) == round(c3.m.val_SI * c3.fluid.val["water"], 2) +assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 2) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 2) -# replace initial draft of constraints -> maybe add saturation constraint as parameter -absorber.get_mandatory_constraints = absorber.new_constraints +c4.set_attr(h=None) +pu.set_attr(eta_s=0.8) # now there is two additional equations in case both fluids at outlet 1 are variable # and one additional equation if none of the fluids is variale # fluid composition as function of saturation temperature and pressure @@ -412,14 +391,43 @@ def new_constraints(self): c3.h.val0 = c3.h.val_SI nw.solve("design") # check results, should be identical to previous ones -assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) -assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) +assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 2) == round(c3.m.val_SI * c3.fluid.val["water"], 2) +assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 2) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 2) # fix the either the water or the libr mass fraction and the temperature or pressure # gives us the other respective value, pressure in this case -c3.set_attr(fluid={"INCOMP::LiBr": 0.45}, T=295, p=None) +c3.set_attr(fluid={"INCOMP::LiBr": 0.50}, T=295, p=None) nw.solve("design") assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) print(f"T3: {c3.T.val_SI}, p3: {c3.p.val_SI}, fluid3: {c3.fluid.val}") + +nw.print_results() + + +he = HeatExchanger("internal heat exchanger") + +nw.del_conns(c4, c7) + +c4 = Connection(pu, "out1", he, "in2", label="4") +c5 = Connection(he, "out2", rich, "in1", label="5") + +c6 = Connection(poor, "out1", he, "in1", label="6") +c7 = Connection(he, "out1", va, "in1", label="7") + +nw.add_conns(c4, c5, c6, c7) + +c4.set_attr(p=p_cond, h0=37400) +c5.set_attr(h0=40000, p=p_cond) +c6.set_attr(h=h_poor, p=p_cond) +c7.set_attr(p=p_cond, h=100000) +c2.set_attr(h=None) + +nw.solve("design") + +he.set_attr(ttd_l=5) +c7.set_attr(h=None) +nw.solve("design") + +nw.print_results() From 08df71d5e5589b4199d82ded4eabd5466c3ef325 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Tue, 31 Oct 2023 12:58:59 +0100 Subject: [PATCH 10/22] Fine tune temperature delta for derivatives --- src/tespy/connections/connection.py | 7 ++- src/tespy/networks/network.py | 2 +- src/tespy/tools/fluid_properties/helpers.py | 6 ++- src/tespy/tools/fluid_properties/mixtures.py | 47 +++++++++----------- 4 files changed, 31 insertions(+), 31 deletions(-) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index eb9d8db19..5842aeeb8 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -910,8 +910,7 @@ def calc_results(self): number_fluids = get_number_of_fluids(self.fluid_data) if number_fluids > 1: h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule, solvent=self.solvent) - print(h_from_T, self.h.val_SI) - if abs(h_from_T - self.h.val_SI) > 1e-1: + if abs(h_from_T - self.h.val_SI) > ERR ** 0.5: self.T.val_SI = np.nan self.vol.val_SI = np.nan self.v.val_SI = np.nan @@ -1029,8 +1028,8 @@ def check_temperature_bounds(self): Tmax = min( [w._T_max for f, w in self.fluid.wrapper.items() if self.fluid.val[f] > ERR] ) * 0.99 - hmin = h_mix_pT(self.p.val_SI, Tmin, self.fluid_data, self.mixing_rule) - hmax = h_mix_pT(self.p.val_SI, Tmax, self.fluid_data, self.mixing_rule) + hmin = h_mix_pT(self.p.val_SI, Tmin, self.fluid_data, self.mixing_rule, solvent=self.solvent) + hmax = h_mix_pT(self.p.val_SI, Tmax, self.fluid_data, self.mixing_rule, solvent=self.solvent) if self.h.val_SI < hmin: self.h.val_SI = hmin diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index bf6f66521..c992f7143 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -1776,7 +1776,7 @@ def init_precalc_properties(self, c): if c.T.is_set: try: - c.h.val_SI = fp.h_mix_pT(c.p.val_SI, c.T.val_SI, c.fluid_data, c.mixing_rule) + c.h.val_SI = fp.h_mix_pT(c.p.val_SI, c.T.val_SI, c.fluid_data, c.mixing_rule, solvent=c.solvent) except ValueError: pass diff --git a/src/tespy/tools/fluid_properties/helpers.py b/src/tespy/tools/fluid_properties/helpers.py index 3b6f34e55..58c8b5296 100644 --- a/src/tespy/tools/fluid_properties/helpers.py +++ b/src/tespy/tools/fluid_properties/helpers.py @@ -70,9 +70,13 @@ def inverse_temperature_mixture(p=None, target_value=None, fluid_data=None, T0=N if T0 is None: T0 = (valmin + valmax) / 2.0 T0 = 320 + if "solvent" in kwargs: + delta = 1e-5 + else: + delta = 1e-2 function_kwargs = { "p": p, "fluid_data": fluid_data, "T": T0, - "function": f, "parameter": "T" , "delta": 0.01 + "function": f, "parameter": "T" , "delta": delta } function_kwargs.update(**kwargs) return newton_with_kwargs( diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index c201f94f7..f018fae0c 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -35,37 +35,37 @@ def h_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): def xsat_pT_incomp_solution(p=None, T=None, fluid_data=None, **kwargs): - X_min = 0 - X_max = 0.75 - X = .3 + x_min = 0 + x_max = 0.75 + x = kwargs["x0"] d = 1e-5 iter = 0 solvent = kwargs["solvent"] while True: - res = fluid_data[solvent]["wrapper"].psat_Tx(T, X) - p - upper = fluid_data[solvent]["wrapper"].psat_Tx(T, X + d) - lower = fluid_data[solvent]["wrapper"].psat_Tx(T, X - d) - fluid_data[solvent]["wrapper"].AS.set_mass_fractions([X]) + res = fluid_data[solvent]["wrapper"].psat_Tx(T, x) - p + upper = fluid_data[solvent]["wrapper"].psat_Tx(T, x + d) + lower = fluid_data[solvent]["wrapper"].psat_Tx(T, x - d) + fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x]) deriv = (upper - lower) / (2 * d) - X -= res / deriv + x -= res / deriv if abs(res) < 1e-6: - if X < X_min: - return X_min - elif X > X_max: - return X_max + if x < x_min: + return x_min + elif x > x_max: + return x_max break - if X >= X_max: - X = X_max - d - elif X <= X_min: - X = X_min + d + if x >= x_max: + x = x_max - d + elif x <= x_min: + x = x_min + d iter += 1 if iter > 10: break - return X + return x def h_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): @@ -123,17 +123,16 @@ def h_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): solvent = kwargs["solvent"] if p < fluid_data[solvent]["wrapper"].p_sat(T): x_old = fluid_data[solvent]["mass_fraction"] + kwargs["x0"] = x_old x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs) + fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_new]) x_water_gas = x_new - x_old h = fluid_data[solvent]["wrapper"].h_pT(p + 1e-6, T) * (1 - x_water_gas) h += fluid_data["water"]["wrapper"].h_QT(1, T) * x_water_gas fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old]) return h else: - try: - return fluid_data[solvent]["wrapper"].h_pT(p, T) - except ValueError: - return fluid_data[solvent]["wrapper"].h_pT(p + 1e-6, T) + return fluid_data[solvent]["wrapper"].h_pT(p, T) def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): @@ -186,6 +185,7 @@ def s_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): solvent = kwargs["solvent"] if p < fluid_data[solvent]["wrapper"].p_sat(T): x_old = fluid_data[solvent]["mass_fraction"] + kwargs["x0"] = x_old x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs) x_water_gas = x_new - x_old s = fluid_data[solvent]["wrapper"].s_pT(p + 1e-6, T) * (1 - x_water_gas) @@ -193,10 +193,7 @@ def s_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old]) return s else: - try: - return fluid_data[solvent]["wrapper"].s_pT(p, T) - except ValueError: - return fluid_data[solvent]["wrapper"].s_pT(p + 1e-6, T) + return fluid_data[solvent]["wrapper"].s_pT(p, T) def v_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): From f697d4c497a8372d8de3cf68f7391efb059efd8d Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Tue, 31 Oct 2023 12:59:19 +0100 Subject: [PATCH 11/22] Fix fluid mass fractions for numerical derivatives towards fluid composition --- src/tespy/tools/helpers.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/tespy/tools/helpers.py b/src/tespy/tools/helpers.py index e177474ec..019f306df 100644 --- a/src/tespy/tools/helpers.py +++ b/src/tespy/tools/helpers.py @@ -444,13 +444,16 @@ def _numeric_deriv(obj, func, dx, conn=None, **kwargs): elif dx in conn.fluid.is_var: d = 1e-5 - + variable_sum = sum(conn.fluid.val[f] for f in conn.fluid.is_var if f != dx) + variable_fluids = {f: x for f, x in conn.fluid.val.items() if f in conn.fluid.is_var and f != dx} val = conn.fluid.val[dx] if conn.fluid.val[dx] + d <= 1: conn.fluid.val[dx] += d else: conn.fluid.val[dx] = 1 + for f in variable_fluids: + conn.fluid.val[f] = variable_fluids[f] * (1 - d / variable_sum) conn.build_fluid_data() exp = func(**kwargs) @@ -459,9 +462,15 @@ def _numeric_deriv(obj, func, dx, conn=None, **kwargs): else: conn.fluid.val[dx] = 0 + for f in variable_fluids: + conn.fluid.val[f] = variable_fluids[f] * (1 + d / variable_sum) + conn.build_fluid_data() exp -= func(**kwargs) + for f in variable_fluids: + conn.fluid.val[f] = variable_fluids[f] + conn.fluid.val[dx] = val conn.build_fluid_data() From a68d0ee56cdde096e644cc5a1a08999724495257 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Tue, 31 Oct 2023 12:59:46 +0100 Subject: [PATCH 12/22] Make a system using an absorber only --- test_libr.py | 461 +++++++++++++++++++++++---------------------------- 1 file changed, 212 insertions(+), 249 deletions(-) diff --git a/test_libr.py b/test_libr.py index 5392fe817..ac7de51ec 100644 --- a/test_libr.py +++ b/test_libr.py @@ -4,169 +4,6 @@ from tespy.components.component import Component - -libr = AbstractState("INCOMP", "LiBr") - -libr.set_mass_fractions([0.5]) -libr.update(CP.QT_INPUTS, 0, 300) -libr.p() -libr.set_mass_fractions([0.2]) -libr.update(CP.QT_INPUTS, 0, 300) -libr.p() - -# saturation pressure by temperature and libr mass fraction x - -def psat_TX(T, X): - libr.set_mass_fractions([X]) - libr.update(CP.QT_INPUTS, 0, T) - return libr.p() - - -def Xsat_Tp(T, p): - X_min = 0 - X_max = 0.75 - X = .3 - d = 1e-5 - while True: - res = psat_TX(T, X) - p - deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) - X -= res / deriv - # print(X) - - if abs(res) < 1e-6: - if X < X_min: - return X_min - elif X > X_max: - return X_max - break - - if X >= X_max: - X = X_max - d - elif X <= X_min: - X = X_min + d - return X - - -def Tsat_pX(p, X): - return newton(psat_TX, p, X) - - -def newton(func, p, X): - d = 1e-3 - # T = 300 - T_min = libr.trivial_keyed_output(CP.iT_min) - T_max = libr.trivial_keyed_output(CP.iT_max) - T = (T_min + T_max) / 2 - while True: - res = func(T, X) - p - deriv = (func(T + d, X) - func(T -d, X)) / (2 * d) - T -= res / deriv - - if T >= T_max: - T = T_max - d - elif T <= T_min: - T = T_min + d - - if abs(res) < 1e-6: - return T - - -# x_range = np.linspace(0, 0.75, 76) -# for x in x_range: -# psat_TX(300, x) -# Tsat_pX(1e4, x) - - -def Xsat_Tp(T, p): - X_min = 0 - X_max = 0.75 - X = .3 - d = 1e-5 - while True: - res = psat_TX(T, X) - p - deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) - X -= res / deriv - # print(X) - - if abs(res) < 1e-6: - if X < X_min: - return X_min - elif X > X_max: - return X_max - break - - if X >= X_max: - X = X_max - d - elif X <= X_min: - X = X_min + d - return X - -# X = Xsat_Tp(350, 1e4) -# libr.set_mass_fractions([X]) -# print(X) - -# p_max = psat_TX(400, 0) -# print(p_max, Xsat_Tp(400, p_max)) -# for p in np.linspace(1e4, p_max): -# # print(p) -# print(p, Xsat_Tp(400, p)) - -# libr.update(CP.PT_INPUTS, 1e4, 400) - -# Tsat_pX(1e4, .3) -# libr -# T = libr.s -# libr.update(CP.QT_INPUTS, 0, T) - -# print(libr.trivial_keyed_output(CP.iT_min)) -# print(libr.trivial_keyed_output(CP.iT_freeze)) -# print(libr.trivial_keyed_output(CP.iP)) - -# watercycle -T_evap = 275.15 -p_evap = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_evap, "water") - -T_cond = 309.15 -p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") - -# heat source temperature -T_sol_des_out = 350 - -# heat sink temperature -T_sol_abs_out = 306.15 - -x_rich = Xsat_Tp(T_sol_abs_out, p_evap) - -x_water_rich = 1 - x_rich -T_water_abs_in = T_evap -h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") - -h_sol_abs_out = libr.hmass() -s_sol_abs_out = libr.smass() - -eta_s_pump = 0.9 -libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) -h_pump_out_s = libr.hmass() -h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump - -x_poor = Xsat_Tp(T_sol_des_out, p_cond) -h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) - -delta_T = 0 -T_water_des_out = T_sol_des_out - delta_T -h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") - -x_water_poor = 1 - x_poor - -m_water = 1 -m_rich = m_water * (1 - x_water_poor) / (x_water_rich - x_water_poor) - -m_poor = m_rich - m_water - -delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out -delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out - - from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution @@ -213,7 +50,7 @@ def get_mandatory_constraints(self): "num_eq": 1 }, } - if "LiBr" in self.outl[0].fluid.is_var: + if self.outl[0].solvent in self.outl[0].fluid.is_var: constraints["saturation_constraints_water"] = { "func": self.saturated_solution_water_func, "deriv": self.saturated_solution_water_deriv, @@ -235,18 +72,22 @@ def mass_flow_deriv(self, k): def fluid_func(self): return ( - self.inl[0].m.val_SI - + self.inl[1].m.val_SI * self.inl[1].fluid.val[self.h2o] - - self.outl[0].m.val_SI * self.outl[0].fluid.val[self.h2o] + self.inl[1].m.val_SI * self.inl[1].fluid.val["LiBr"] + - self.outl[0].m.val_SI * self.outl[0].fluid.val["LiBr"] ) def fluid_deriv(self, increment_filter, k): - if self.inl[0].m.is_var: - self.jacobian[k, self.inl[0].m.J_col] = 1 - if self.inl[1].m.is_var: - self.jacobian[k, self.inl[1].m.J_col] = self.inl[1].fluid.val[self.h2o] - if self.outl[0].m.is_var: - self.jacobian[k, self.outl[0].m.J_col] = -self.outl[0].fluid.val[self.h2o] + outl = self.outl[0] + inl = self.inl[1] + + if inl.m.is_var: + self.jacobian[k, inl.m.J_col] = inl.fluid.val[inl.solvent] + if inl.solvent in inl.fluid.is_var: + self.jacobian[k, inl.fluid.J_col[inl.solvent]] = inl.m.val_SI + if outl.m.is_var: + self.jacobian[k, outl.m.J_col] = -outl.fluid.val[outl.solvent] + if outl.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -outl.m.val_SI def pressure_equality_func(self): residual = [] @@ -287,7 +128,7 @@ def saturated_solution_libr_func(self): outl = self.outl[0] x_previous = outl.fluid.val[outl.solvent] T = outl.calc_T() - x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=outl.solvent) + x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=outl.solvent, x0=x_previous) outl.fluid_data[outl.solvent]["wrapper"].AS.set_mass_fractions([x_previous]) return x_libr - outl.fluid.val[outl.solvent] @@ -300,7 +141,7 @@ def saturated_solution_libr_deriv(self, increment_filter, k): deriv = self.numeric_deriv(self.saturated_solution_libr_func, "h", outl) self.jacobian[k, outl.h.J_col] = deriv if outl.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -1 + self.jacobian[k, outl.fluid.J_col[outl.solvent]] = self.numeric_deriv(self.saturated_solution_libr_func, outl.solvent, outl) @staticmethod def is_branch_source(): @@ -332,102 +173,224 @@ def propagate_wrapper_to_target(self, branch): outconn.target.propagate_wrapper_to_target(branch) -from tespy.components import Source, Sink, Pump, Valve, HeatExchanger -from tespy.networks import Network -from tespy.connections import Connection +if __name__ == "__main__": + libr = AbstractState("INCOMP", "LiBr") + + libr.set_mass_fractions([0.5]) + libr.update(CP.QT_INPUTS, 0, 300) + libr.p() + libr.set_mass_fractions([0.2]) + libr.update(CP.QT_INPUTS, 0, 300) + libr.p() + + # saturation pressure by temperature and libr mass fraction x + + def psat_TX(T, X): + libr.set_mass_fractions([X]) + libr.update(CP.QT_INPUTS, 0, T) + return libr.p() + + + def Xsat_Tp(T, p): + X_min = 0 + X_max = 0.75 + X = .3 + d = 1e-5 + while True: + res = psat_TX(T, X) - p + deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) + X -= res / deriv + # print(X) + + if abs(res) < 1e-6: + if X < X_min: + return X_min + elif X > X_max: + return X_max + break + + if X >= X_max: + X = X_max - d + elif X <= X_min: + X = X_min + d + + psat_TX(T, X) + return X + + + def Tsat_pX(p, X): + return newton(psat_TX, p, X) + + + def newton(func, p, X): + d = 1e-3 + # T = 300 + T_min = libr.trivial_keyed_output(CP.iT_min) + T_max = libr.trivial_keyed_output(CP.iT_max) + T = (T_min + T_max) / 2 + while True: + res = func(T, X) - p + deriv = (func(T + d, X) - func(T -d, X)) / (2 * d) + T -= res / deriv + + if T >= T_max: + T = T_max - d + elif T <= T_min: + T = T_min + d + + if abs(res) < 1e-6: + return T + + + # X = Xsat_Tp(350, 1e4) + # libr.set_mass_fractions([X]) + # print(X) + + # p_max = psat_TX(400, 0) + # print(p_max, Xsat_Tp(400, p_max)) + # for p in np.linspace(1e4, p_max): + # # print(p) + # print(p, Xsat_Tp(400, p)) + + # libr.update(CP.PT_INPUTS, 1e4, 400) + + # Tsat_pX(1e4, .3) + # libr + # T = libr.s + # libr.update(CP.QT_INPUTS, 0, T) + + # print(libr.trivial_keyed_output(CP.iT_min)) + # print(libr.trivial_keyed_output(CP.iT_freeze)) + # print(libr.trivial_keyed_output(CP.iP)) + + +# Absorber +# +# in1: water +# in2: wasser/LiBr Lösung mit geringem Wasseranteil bzw. Mischung aus Lösung und Wasserdampf +# out1: gesättigte LiBr Lösung +# m_in1 + m_in2 = m_out1 +# m_in1 * x_water,in1 + m_in2 * x_water,in2 = m_out1 * x_water,out1 +# x_libr,out1 = x(p_out1, T_out1, voll gesättigt) p_sat(x, T_sat) berechnen +# x_water,out1 = 1 - x_libr,out1 +# p_in1 = p_out1 +# p_in2 = p_out2 +# + # watercycle + T_evap = 275.15 + p_evap = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_evap, "water") + + T_cond = 309.15 + p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") + + # heat source temperature + T_sol_des_out = 350 + + # heat sink temperature + T_sol_abs_out = 306.15 + + x_rich = Xsat_Tp(T_sol_abs_out, p_evap) + + x_water_rich = 1 - x_rich + T_water_abs_in = T_evap + h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") + + h_sol_abs_out = libr.hmass() + s_sol_abs_out = libr.smass() + + eta_s_pump = 0.9 + libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) + h_pump_out_s = libr.hmass() + h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump + x_poor = Xsat_Tp(T_sol_des_out, p_cond) + h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) -nw = Network() + delta_T = 0 + T_water_des_out = T_sol_des_out - delta_T + h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") -water = Source("water source") -rich = Sink("rich solution") -poor = Source("poor solution") + x_water_poor = 1 - x_poor -pu = Pump("pump") -va = Valve("valve") + m_water = 1 + m_rich = m_water * (1 - x_water_poor) / (x_water_rich - x_water_poor) -absorber = Absorber("absorber") + m_poor = m_rich - m_water -c1 = Connection(water, "out1", absorber, "in1", label="1") -c2 = Connection(va, "out1", absorber, "in2", label="2") -c3 = Connection(absorber, "out1", pu, "in1", label="3") -c4 = Connection(pu, "out1", rich, "in1", label="4") -c7 = Connection(poor, "out1", va, "in1", label="7") + delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out + delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out -nw.add_conns(c1, c2, c3, c4, c7) + from tespy.components import Source, Sink + from tespy.networks import Network + from tespy.connections import Connection -# pu.set_attr(eta_s=0.8) -c1.set_attr(fluid={"water": 1}, p=p_evap, m=1, x=1) -c2.set_attr(fluid={"water": x_water_poor, "INCOMP::LiBr": 1 - x_water_poor}, h=h_poor, mixing_rule="incomp-solution", solvent="LiBr") -c3.set_attr(h=h_sol_abs_out, p0=p_evap) -print(h_sol_abs_out) -c4.set_attr(p=p_cond, h=37410) -c7.set_attr(p=p_cond) -nw.set_attr() -nw.solve("design") + nw = Network() -# how does the reference temperature affect the energy balance -# p_ref = 1e5 -# for T_ref in range(274, 400): -# h_water_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, "water") -# h_poor_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_poor}]") -# h_rich_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{x_rich}]") -# print((c1.h.val_SI - h_water_ref) * c1.m.val_SI + (c2.h.val_SI - h_poor_ref) * c2.m.val_SI - (c3.h.val_SI - h_rich_ref) * c3.m.val_SI) + water = Source("water source") + rich = Sink("rich solution") + poor = Source("poor solution") -# some checks -nw.print_results() -assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 2) == round(c3.m.val_SI * c3.fluid.val["water"], 2) -assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 2) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 2) + absorber = Absorber("absorber") -c4.set_attr(h=None) -pu.set_attr(eta_s=0.8) -# now there is two additional equations in case both fluids at outlet 1 are variable -# and one additional equation if none of the fluids is variale -# fluid composition as function of saturation temperature and pressure -c3.fluid.is_set = set() -c1.set_attr(p=None) -c3.set_attr(h=None, T=306.15, p=p_evap) -c3.h.val0 = c3.h.val_SI -nw.solve("design") -# check results, should be identical to previous ones -assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 2) == round(c3.m.val_SI * c3.fluid.val["water"], 2) -assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 2) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 2) + c1 = Connection(water, "out1", absorber, "in1", label="1") + c2 = Connection(poor, "out1", absorber, "in2", label="2") + c3 = Connection(absorber, "out1", rich, "in1", label="3") -# fix the either the water or the libr mass fraction and the temperature or pressure -# gives us the other respective value, pressure in this case -c3.set_attr(fluid={"INCOMP::LiBr": 0.50}, T=295, p=None) -nw.solve("design") -assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) -assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) + nw.add_conns(c1, c2, c3) -print(f"T3: {c3.T.val_SI}, p3: {c3.p.val_SI}, fluid3: {c3.fluid.val}") + c1.set_attr(fluid={"water": 1}, p=0.01e5, m=1, x=1) + c2.set_attr(fluid={"water": x_water_poor, "INCOMP::LiBr": 1 - x_water_poor}, h=h_poor, mixing_rule="incomp-solution", solvent="LiBr") + c3.set_attr(h=h_sol_abs_out, p0=0.01e5, m0=m_rich) -nw.print_results() + nw.solve("design") + c2.set_attr(m=10) + c3.set_attr(h=None) -he = HeatExchanger("internal heat exchanger") + nw.solve("design") + nw.print_results() + print(c3.fluid.val) + print(c2.fluid.val) -nw.del_conns(c4, c7) + print(c1.m.val_SI * c1.h.val_SI + c2.m.val_SI * c2.h.val_SI - c3.m.val_SI * c3.h.val_SI) -c4 = Connection(pu, "out1", he, "in2", label="4") -c5 = Connection(he, "out2", rich, "in1", label="5") + # how does the reference temperature affect the energy balance + # eb_prev = 0 + # p_ref = 1e5 + # for T_ref in range(274, 350): + # h_water_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, "water") + # h_poor_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{c2.fluid.val['LiBr']}]") + # h_rich_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{c3.fluid.val['LiBr']}]") + # (CP.CoolProp.PropsSI("Cpmass", "P", p_ref, "T", T_ref, "water")) + # eb = (c1.h.val_SI - h_water_ref) * c1.m.val_SI + (c2.h.val_SI - h_poor_ref) * c2.m.val_SI - (c3.h.val_SI - h_rich_ref) * c3.m.val_SI + # eb_prev = eb -c6 = Connection(poor, "out1", he, "in1", label="6") -c7 = Connection(he, "out1", va, "in1", label="7") + # s -nw.add_conns(c4, c5, c6, c7) + # some checks + # nw.print_results() + print(m_rich, c3.m.val_SI) + assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) + assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) -c4.set_attr(p=p_cond, h0=37400) -c5.set_attr(h0=40000, p=p_cond) -c6.set_attr(h=h_poor, p=p_cond) -c7.set_attr(p=p_cond, h=100000) -c2.set_attr(h=None) + # use temperature instead of enthalpy + c2.set_attr(m=None) + c3.set_attr(T=306.15) + nw.solve("design") + nw.print_results() + # check results, should be identical to previous ones + assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) + assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) -nw.solve("design") + # fix the either the water or the libr mass fraction and the temperature or pressure + # gives us the other respective value, pressure in this case + c1.set_attr(p=None) + c3.set_attr(fluid={"INCOMP::LiBr": 0.50}, T=305) + nw.solve("design") + assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) + assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) -he.set_attr(ttd_l=5) -c7.set_attr(h=None) -nw.solve("design") + print(f"T3: {c3.T.val_SI}, p3: {c3.p.val_SI}, fluid3: {c3.fluid.val}") -nw.print_results() + nw.print_results() From ebc75c62c9f3dea4f2ca1725f30d3c4931b84595 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Tue, 31 Oct 2023 13:00:07 +0100 Subject: [PATCH 13/22] Rename file --- test_libr.py => test_absorber.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test_libr.py => test_absorber.py (100%) diff --git a/test_libr.py b/test_absorber.py similarity index 100% rename from test_libr.py rename to test_absorber.py From 06e2438c341dbdebba692d1839ef2da39ec65dec Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Tue, 31 Oct 2023 13:07:12 +0100 Subject: [PATCH 14/22] Separate reference calculation from absorber implementation --- sorption_reference.py | 116 ++++++++++++++++++++++++++++++++ test_absorber.py | 151 +----------------------------------------- 2 files changed, 118 insertions(+), 149 deletions(-) create mode 100644 sorption_reference.py diff --git a/sorption_reference.py b/sorption_reference.py new file mode 100644 index 000000000..616e084d0 --- /dev/null +++ b/sorption_reference.py @@ -0,0 +1,116 @@ + +from CoolProp.CoolProp import AbstractState +import CoolProp as CP +import numpy as np + + +libr = AbstractState("INCOMP", "LiBr") + +libr.set_mass_fractions([0.5]) +libr.update(CP.QT_INPUTS, 0, 300) +libr.p() +libr.set_mass_fractions([0.2]) +libr.update(CP.QT_INPUTS, 0, 300) +libr.p() + +# saturation pressure by temperature and libr mass fraction x + +def psat_TX(T, X): + libr.set_mass_fractions([X]) + libr.update(CP.QT_INPUTS, 0, T) + return libr.p() + + +def Xsat_Tp(T, p): + X_min = 0 + X_max = 0.75 + X = .3 + d = 1e-5 + while True: + res = psat_TX(T, X) - p + deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) + X -= res / deriv + # print(X) + + if abs(res) < 1e-6: + if X < X_min: + return X_min + elif X > X_max: + return X_max + break + + if X >= X_max: + X = X_max - d + elif X <= X_min: + X = X_min + d + + psat_TX(T, X) + return X + + +def Tsat_pX(p, X): + return newton(psat_TX, p, X) + + +def newton(func, p, X): + d = 1e-3 + # T = 300 + T_min = libr.trivial_keyed_output(CP.iT_min) + T_max = libr.trivial_keyed_output(CP.iT_max) + T = (T_min + T_max) / 2 + while True: + res = func(T, X) - p + deriv = (func(T + d, X) - func(T -d, X)) / (2 * d) + T -= res / deriv + + if T >= T_max: + T = T_max - d + elif T <= T_min: + T = T_min + d + + if abs(res) < 1e-6: + return T + +# watercycle +T_evap = 275.15 +p_evap = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_evap, "water") + +T_cond = 309.15 +p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") + +# heat source temperature +T_sol_des_out = 350 + +# heat sink temperature +T_sol_abs_out = 306.15 + +x_rich = Xsat_Tp(T_sol_abs_out, p_evap) + +x_water_rich = 1 - x_rich +T_water_abs_in = T_evap +h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") + +h_sol_abs_out = libr.hmass() +s_sol_abs_out = libr.smass() + +eta_s_pump = 0.9 +libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) +h_pump_out_s = libr.hmass() +h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump + +x_poor = Xsat_Tp(T_sol_des_out, p_cond) +h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) + +delta_T = 0 +T_water_des_out = T_sol_des_out - delta_T +h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") + +x_water_poor = 1 - x_poor + +m_water = 1 +m_rich = m_water * (1 - x_water_poor) / (x_water_rich - x_water_poor) + +m_poor = m_rich - m_water + +delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out +delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out diff --git a/test_absorber.py b/test_absorber.py index ac7de51ec..03fb18007 100644 --- a/test_absorber.py +++ b/test_absorber.py @@ -1,10 +1,9 @@ -from CoolProp.CoolProp import AbstractState -import CoolProp as CP -import numpy as np from tespy.components.component import Component from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution +# this runs the references calculations and imports all variables from there +from sorption_reference import * class Absorber(Component): @@ -122,7 +121,6 @@ def saturated_solution_water_deriv(self, increment_filter, k): self.jacobian[k, outl.fluid.J_col[self.h2o]] = -1 if outl.solvent in outl.fluid.is_var: self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -1 - # pass def saturated_solution_libr_func(self): outl = self.outl[0] @@ -174,151 +172,6 @@ def propagate_wrapper_to_target(self, branch): if __name__ == "__main__": - libr = AbstractState("INCOMP", "LiBr") - - libr.set_mass_fractions([0.5]) - libr.update(CP.QT_INPUTS, 0, 300) - libr.p() - libr.set_mass_fractions([0.2]) - libr.update(CP.QT_INPUTS, 0, 300) - libr.p() - - # saturation pressure by temperature and libr mass fraction x - - def psat_TX(T, X): - libr.set_mass_fractions([X]) - libr.update(CP.QT_INPUTS, 0, T) - return libr.p() - - - def Xsat_Tp(T, p): - X_min = 0 - X_max = 0.75 - X = .3 - d = 1e-5 - while True: - res = psat_TX(T, X) - p - deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) - X -= res / deriv - # print(X) - - if abs(res) < 1e-6: - if X < X_min: - return X_min - elif X > X_max: - return X_max - break - - if X >= X_max: - X = X_max - d - elif X <= X_min: - X = X_min + d - - psat_TX(T, X) - return X - - - def Tsat_pX(p, X): - return newton(psat_TX, p, X) - - - def newton(func, p, X): - d = 1e-3 - # T = 300 - T_min = libr.trivial_keyed_output(CP.iT_min) - T_max = libr.trivial_keyed_output(CP.iT_max) - T = (T_min + T_max) / 2 - while True: - res = func(T, X) - p - deriv = (func(T + d, X) - func(T -d, X)) / (2 * d) - T -= res / deriv - - if T >= T_max: - T = T_max - d - elif T <= T_min: - T = T_min + d - - if abs(res) < 1e-6: - return T - - - # X = Xsat_Tp(350, 1e4) - # libr.set_mass_fractions([X]) - # print(X) - - # p_max = psat_TX(400, 0) - # print(p_max, Xsat_Tp(400, p_max)) - # for p in np.linspace(1e4, p_max): - # # print(p) - # print(p, Xsat_Tp(400, p)) - - # libr.update(CP.PT_INPUTS, 1e4, 400) - - # Tsat_pX(1e4, .3) - # libr - # T = libr.s - # libr.update(CP.QT_INPUTS, 0, T) - - # print(libr.trivial_keyed_output(CP.iT_min)) - # print(libr.trivial_keyed_output(CP.iT_freeze)) - # print(libr.trivial_keyed_output(CP.iP)) - - -# Absorber -# -# in1: water -# in2: wasser/LiBr Lösung mit geringem Wasseranteil bzw. Mischung aus Lösung und Wasserdampf -# out1: gesättigte LiBr Lösung -# m_in1 + m_in2 = m_out1 -# m_in1 * x_water,in1 + m_in2 * x_water,in2 = m_out1 * x_water,out1 -# x_libr,out1 = x(p_out1, T_out1, voll gesättigt) p_sat(x, T_sat) berechnen -# x_water,out1 = 1 - x_libr,out1 -# p_in1 = p_out1 -# p_in2 = p_out2 -# - # watercycle - T_evap = 275.15 - p_evap = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_evap, "water") - - T_cond = 309.15 - p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") - - # heat source temperature - T_sol_des_out = 350 - - # heat sink temperature - T_sol_abs_out = 306.15 - - x_rich = Xsat_Tp(T_sol_abs_out, p_evap) - - x_water_rich = 1 - x_rich - T_water_abs_in = T_evap - h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") - - h_sol_abs_out = libr.hmass() - s_sol_abs_out = libr.smass() - - eta_s_pump = 0.9 - libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) - h_pump_out_s = libr.hmass() - h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump - - x_poor = Xsat_Tp(T_sol_des_out, p_cond) - h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) - - delta_T = 0 - T_water_des_out = T_sol_des_out - delta_T - h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") - - x_water_poor = 1 - x_poor - - m_water = 1 - m_rich = m_water * (1 - x_water_poor) / (x_water_rich - x_water_poor) - - m_poor = m_rich - m_water - - delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out - delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out from tespy.components import Source, Sink from tespy.networks import Network From 06b8061354598b8e4a69e6a7e103c641c0919fd2 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Thu, 2 Nov 2023 19:42:09 +0100 Subject: [PATCH 15/22] Add a method to extract the fluid wrapper sources from network's components --- src/tespy/components/basics/cycle_closer.py | 4 ++++ src/tespy/components/basics/source.py | 4 ++++ src/tespy/components/component.py | 4 ++++ src/tespy/components/reactors/fuel_cell.py | 6 ++++-- src/tespy/components/reactors/water_electrolyzer.py | 4 ++++ src/tespy/networks/network.py | 4 +--- 6 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/tespy/components/basics/cycle_closer.py b/src/tespy/components/basics/cycle_closer.py index 23aba26da..754b2c75f 100644 --- a/src/tespy/components/basics/cycle_closer.py +++ b/src/tespy/components/basics/cycle_closer.py @@ -131,6 +131,10 @@ def outlets(): def is_branch_source(): return True + @staticmethod + def is_wrapper_branch_source(): + return True + def start_branch(self): outconn = self.outl[0] branch = { diff --git a/src/tespy/components/basics/source.py b/src/tespy/components/basics/source.py index 35719fa34..8f1598c0e 100644 --- a/src/tespy/components/basics/source.py +++ b/src/tespy/components/basics/source.py @@ -69,6 +69,10 @@ def outlets(): def is_branch_source(): return True + @staticmethod + def is_wrapper_branch_source(): + return True + def start_branch(self): outconn = self.outl[0] branch = { diff --git a/src/tespy/components/component.py b/src/tespy/components/component.py index 8c8dec496..5a90ef1b4 100644 --- a/src/tespy/components/component.py +++ b/src/tespy/components/component.py @@ -309,6 +309,10 @@ def _serializable(): def is_branch_source(): return False + @staticmethod + def is_wrapper_branch_source(): + return False + def propagate_to_target(self, branch): inconn = branch["connections"][-1] conn_idx = self.inl.index(inconn) diff --git a/src/tespy/components/reactors/fuel_cell.py b/src/tespy/components/reactors/fuel_cell.py index 088d09757..3c000c97e 100644 --- a/src/tespy/components/reactors/fuel_cell.py +++ b/src/tespy/components/reactors/fuel_cell.py @@ -708,12 +708,14 @@ def calc_P(self): ) return val - - @staticmethod def is_branch_source(): return True + @staticmethod + def is_wrapper_branch_source(): + return True + def start_branch(self): outconn = self.outl[1] if "H2O" not in outconn.fluid.val: diff --git a/src/tespy/components/reactors/water_electrolyzer.py b/src/tespy/components/reactors/water_electrolyzer.py index 9e8e959d8..fec371aeb 100644 --- a/src/tespy/components/reactors/water_electrolyzer.py +++ b/src/tespy/components/reactors/water_electrolyzer.py @@ -1037,6 +1037,10 @@ def bus_deriv(self, bus): def is_branch_source(): return True + @staticmethod + def is_wrapper_branch_source(): + return True + def start_branch(self): branches = {} for outconn in self.outl[1:]: diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index c992f7143..ede0271df 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -724,9 +724,7 @@ def create_massflow_and_fluid_branches(self): def create_fluid_wrapper_branches(self): self.fluid_wrapper_branches = {} - mask = self.comps["comp_type"].isin( - ["Source", "CycleCloser", "WaterElectrolyzer", "FuelCell"] - ) + mask = self.comps["object"].apply(lambda c: c.is_wrapper_branch_source()) start_components = self.comps["object"].loc[mask] for start in start_components: From ad92307026d68091a7b650c983f364a1a241925a Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Thu, 2 Nov 2023 19:42:33 +0100 Subject: [PATCH 16/22] Log deviation in enthalpy in connection property postprocessing --- src/tespy/connections/connection.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index 5842aeeb8..2e0a06590 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -918,7 +918,8 @@ def calc_results(self): msg = ( "Could not find a feasible value for mixture temperature at " f"connection {self.label}. The values for temperature, " - "specific volume, volumetric flow and entropy are set to nan." + "specific volume, volumetric flow and entropy are set to nan. " + f"The deviation is {h_from_T - self.h.val_SI} J/kg." ) logger.error(msg) _converged = False From 4538225100705549855d822bff592759f9243cae Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Thu, 2 Nov 2023 19:43:36 +0100 Subject: [PATCH 17/22] Move import and fix minor bug --- test_absorber.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test_absorber.py b/test_absorber.py index 03fb18007..587664cae 100644 --- a/test_absorber.py +++ b/test_absorber.py @@ -2,8 +2,6 @@ from tespy.components.component import Component from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution -# this runs the references calculations and imports all variables from there -from sorption_reference import * class Absorber(Component): @@ -168,6 +166,7 @@ def propagate_wrapper_to_target(self, branch): outconn = self.outl[0] branch["connections"] += [outconn] + branch["components"] += [self] outconn.target.propagate_wrapper_to_target(branch) @@ -177,6 +176,8 @@ def propagate_wrapper_to_target(self, branch): from tespy.networks import Network from tespy.connections import Connection + # this runs the references calculations and imports all variables from there + from sorption_reference import * nw = Network() @@ -197,14 +198,14 @@ def propagate_wrapper_to_target(self, branch): c3.set_attr(h=h_sol_abs_out, p0=0.01e5, m0=m_rich) nw.solve("design") + nw.print_results() + print(c3.fluid.val) c2.set_attr(m=10) c3.set_attr(h=None) nw.solve("design") nw.print_results() - print(c3.fluid.val) - print(c2.fluid.val) print(c1.m.val_SI * c1.h.val_SI + c2.m.val_SI * c2.h.val_SI - c3.m.val_SI * c3.h.val_SI) From 7e356e4edc589a23c22087f8648142f33ee9d48d Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Thu, 2 Nov 2023 19:44:48 +0100 Subject: [PATCH 18/22] Add desorber model draft --- test_desorber.py | 255 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 255 insertions(+) create mode 100644 test_desorber.py diff --git a/test_desorber.py b/test_desorber.py new file mode 100644 index 000000000..904a092b6 --- /dev/null +++ b/test_desorber.py @@ -0,0 +1,255 @@ + +from tespy.components.component import Component + +from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution + + +class Desorber(Component): + + def __init__(self, label, **kwargs): + super().__init__(label, **kwargs) + + @staticmethod + def inlets(): + return ["in1"] + + @staticmethod + def outlets(): + return ["out1", "out2"] + + def preprocess(self, num_vars): + self.h2o = "water" + super().preprocess(num_vars) + + def get_parameters(self): + return {} + + def get_mandatory_constraints(self): + constraints = { + 'mass_flow_constraints': { + 'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv, + 'constant_deriv': True,# 'latex': self.mass_flow_func_doc, + 'num_eq': 1}, + 'fluid_constraints': { + 'func': self.fluid_func, 'deriv': self.fluid_deriv, + 'constant_deriv': False,# 'latex': self.fluid_func_doc, + 'num_eq': 1}, + 'pressure_constraints': { + 'func': self.pressure_equality_func, + 'deriv': self.pressure_equality_deriv, + 'constant_deriv': True, + 'latex': self.pressure_equality_func_doc, + 'num_eq': 2}, + "saturation_constraints_libr": { + "func": self.saturated_solution_libr_func, + "deriv": self.saturated_solution_libr_deriv, + "constant_deriv": False, + "num_eq": 1 + }, + } + if self.outl[1].solvent in self.outl[1].fluid.is_var: + constraints["saturation_constraints_water"] = { + "func": self.saturated_solution_water_func, + "deriv": self.saturated_solution_water_deriv, + "constant_deriv": False, + "num_eq": 1 + } + return constraints + + def mass_flow_func(self): + return self.inl[0].m.val_SI - self.outl[0].m.val_SI - self.outl[1].m.val_SI + + def mass_flow_deriv(self, k): + inl = self.inl[0] + if inl.m.is_var: + self.jacobian[k, inl.m.J_col] = 1 + + for outl in self.outl: + if outl.m.is_var: + self.jacobian[k, outl.m.J_col] = -1 + + def fluid_func(self): + inl = self.inl[0] + outl = self.outl[1] + return ( + inl.m.val_SI * inl.fluid.val[inl.solvent] + - outl.m.val_SI * outl.fluid.val[outl.solvent] + ) + + def fluid_deriv(self, increment_filter, k): + outl = self.outl[1] + inl = self.inl[0] + + if inl.m.is_var: + self.jacobian[k, inl.m.J_col] = inl.fluid.val[inl.solvent] + if inl.solvent in inl.fluid.is_var: + self.jacobian[k, inl.fluid.J_col[inl.solvent]] = inl.m.val_SI + if outl.m.is_var: + self.jacobian[k, outl.m.J_col] = -outl.fluid.val[outl.solvent] + if outl.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -outl.m.val_SI + + def pressure_equality_func(self): + residual = [] + for c in self.outl : + residual += [c.p.val_SI - self.inl[0].p.val_SI] + return residual + + def pressure_equality_deriv(self, k): + r""" + Calculate partial derivatives for all pressure equations. + + Returns + ------- + deriv : ndarray + Matrix with partial derivatives for the fluid equations. + """ + for c in self.outl: + if c.p.is_var: + self.jacobian[k, c.p.J_col] = 1 + if self.inl[0].p.is_var: + self.jacobian[k, self.inl[0].p.J_col] = -1 + k += 1 + + def saturated_solution_water_func(self): + outl = self.outl[1] + return 1 - outl.fluid.val[outl.solvent] - outl.fluid.val[self.h2o] + + def saturated_solution_water_deriv(self, increment_filter, k): + outl = self.outl[1] + if self.h2o in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[self.h2o]] = -1 + if outl.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -1 + + def saturated_solution_libr_func(self): + outl = self.outl[1] + x_previous = outl.fluid.val[outl.solvent] + T = outl.calc_T() + x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=outl.solvent, x0=x_previous) + outl.fluid_data[outl.solvent]["wrapper"].AS.set_mass_fractions([x_previous]) + return x_libr - outl.fluid.val[outl.solvent] + + def saturated_solution_libr_deriv(self, increment_filter, k): + outl = self.outl[1] + if outl.p.is_var: + deriv = self.numeric_deriv(self.saturated_solution_libr_func, "p", outl) + self.jacobian[k, outl.p.J_col] = deriv + if outl.h.is_var: + deriv = self.numeric_deriv(self.saturated_solution_libr_func, "h", outl) + self.jacobian[k, outl.h.J_col] = deriv + if outl.solvent in outl.fluid.is_var: + self.jacobian[k, outl.fluid.J_col[outl.solvent]] = self.numeric_deriv(self.saturated_solution_libr_func, outl.solvent, outl) + + @staticmethod + def is_branch_source(): + return True + + @staticmethod + def is_wrapper_branch_source(): + return True + + def start_branch(self): + branches = {} + for outconn in self.outl: + branch = { + "connections": [outconn], + "components": [self, outconn.target], + "subbranches": {} + } + outconn.target.propagate_to_target(branch) + + branches[outconn.label] = branch + return branches + + def start_fluid_wrapper_branch(self): + outconn = self.outl[0] + branch = { + "connections": [outconn], + "components": [self] + } + outconn.target.propagate_wrapper_to_target(branch) + + return {outconn.label: branch} + + def propagate_to_target(self, branch): + return + + def propagate_wrapper_to_target(self, branch): + if self in branch["components"]: + return + + outconn = self.outl[1] + branch["connections"] += [outconn] + branch["components"] += [self] + outconn.target.propagate_wrapper_to_target(branch) + + +if __name__ == "__main__": + + from tespy.components import Source, Sink, SimpleHeatExchanger + from tespy.networks import Network + from tespy.connections import Connection + + # this runs the references calculations and imports all variables from there + from sorption_reference import * + + + nw = Network() + + water = Sink("water source") + rich = Source("rich solution") + poor = Sink("poor solution") + + desorber = Desorber("desorber") + valve = SimpleHeatExchanger("valve") + + c5 = Connection(rich, "out1", desorber, "in1", label="5") + c10 = Connection(desorber, "out1", water, "in1", label="10") + c6 = Connection(desorber, "out2", valve, "in1", label="6") + c7 = Connection(valve, "out1", poor, "in1", label="7") + + nw.add_conns(c5, c10, c6, c7) + + c5.set_attr(fluid={"water": x_water_rich, "INCOMP::LiBr": 1 - x_water_rich}, p=p_cond, m=10, h=h_pump_out, mixing_rule="incomp-solution", solvent="LiBr") + c10.set_attr(fluid={"water": 1}, x=1) + c6.set_attr(h=h_poor, p0=0.1e5, m0=m_rich) + c7.set_attr(p=p_cond, h=h_poor * 1.2) + + nw.solve("design") + print(c6.fluid.val) + print(c10.h.val_SI * c10.m.val_SI + c6.h.val_SI * c6.m.val_SI - c5.m.val_SI * c5.h.val_SI) + print(c7.h.val_SI) + print(c7.T.val_SI) + + c10.set_attr(x=None, T=320) + nw.solve("design") + print(c6.fluid.val) + print(c10.h.val_SI * c10.m.val_SI + c6.h.val_SI * c6.m.val_SI - c5.m.val_SI * c5.h.val_SI) + exit() + nw.print_results() + + print(T_sol_des_out) + print(h_poor) + print(c6.h.val_SI) + print(c6.T.val_SI) + print(c10.T.val_SI) + print(c6.fluid.val) + c6.set_attr(m=9) + c6.set_attr(h=None) + + nw.solve("design") + nw.print_results() + print(c6.fluid.val) + + c6.set_attr(m=None, T=350) + nw.solve("design") + nw.print_results() + print(c6.fluid.val) + + c6.set_attr(fluid={"INCOMP::LiBr": 0.7}) + c5.set_attr(p=None) + + nw.solve("design") + nw.print_results() + print(c6.fluid.val) From 46b451e8c25e6288bd40ee1faadcd5d64bb28a3c Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Sun, 5 Nov 2023 10:59:51 +0100 Subject: [PATCH 19/22] Fix fluid wrapper branch definition for desorber --- test_absorber.py | 1 + test_desorber.py | 18 ++++++++++-------- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/test_absorber.py b/test_absorber.py index 587664cae..3425ec943 100644 --- a/test_absorber.py +++ b/test_absorber.py @@ -159,6 +159,7 @@ def propagate_to_target(self, branch): def propagate_wrapper_to_target(self, branch): if branch["connections"][-1] == self.inl[0]: + branch["components"] += [self] return if self in branch["components"]: diff --git a/test_desorber.py b/test_desorber.py index 904a092b6..152f159d3 100644 --- a/test_desorber.py +++ b/test_desorber.py @@ -91,7 +91,7 @@ def fluid_deriv(self, increment_filter, k): def pressure_equality_func(self): residual = [] - for c in self.outl : + for c in self.outl: residual += [c.p.val_SI - self.inl[0].p.val_SI] return residual @@ -163,14 +163,16 @@ def start_branch(self): return branches def start_fluid_wrapper_branch(self): - outconn = self.outl[0] - branch = { - "connections": [outconn], - "components": [self] - } - outconn.target.propagate_wrapper_to_target(branch) + branches = {} + for outconn in self.outl: + branch = { + "connections": [outconn], + "components": [self] + } + outconn.target.propagate_wrapper_to_target(branch) + branches[outconn.label] = branch - return {outconn.label: branch} + return branches def propagate_to_target(self, branch): return From 545d64ebb0f740a72da75f6b066aafb5c5446ac7 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Wed, 14 Feb 2024 08:34:14 +0100 Subject: [PATCH 20/22] Apply some fixes and add volumetric efficiency for pump --- src/tespy/components/turbomachinery/pump.py | 25 ++++ src/tespy/connections/connection.py | 39 +++--- src/tespy/networks/network.py | 9 +- src/tespy/tools/fluid_properties/functions.py | 6 +- src/tespy/tools/fluid_properties/mixtures.py | 21 ++- src/tespy/tools/helpers.py | 25 ++-- test_absorber.py | 123 ++--------------- test_desorber.py | 126 ++---------------- 8 files changed, 116 insertions(+), 258 deletions(-) diff --git a/src/tespy/components/turbomachinery/pump.py b/src/tespy/components/turbomachinery/pump.py index 32d68766e..bd049262f 100644 --- a/src/tespy/components/turbomachinery/pump.py +++ b/src/tespy/components/turbomachinery/pump.py @@ -169,6 +169,11 @@ def get_parameters(self): deriv=self.eta_s_deriv, func=self.eta_s_func, latex=self.eta_s_func_doc), + 'eta_vol': dc_cp( + min_val=0, max_val=1, num_eq=1, + deriv=self.eta_volumetric_deriv, + func=self.eta_volumetric_func, + latex=self.eta_s_func_doc), 'pr': dc_cp( min_val=1, num_eq=1, deriv=self.pr_deriv, @@ -345,6 +350,26 @@ def eta_s_char_deriv(self, increment_filter, k): if self.is_variable(o.h, increment_filter): self.jacobian[k, o.h.J_col] = self.numeric_deriv(f, 'h', o) + def eta_volumetric_func(self): + return ( + self.inl[0].calc_vol() * (self.outl[0].p.val_SI - self.inl[0].p.val_SI) + - (self.outl[0].h.val_SI - self.inl[0].h.val_SI) * self.eta_vol.val + ) + + def eta_volumetric_deriv(self, increment_filter, k): + i = self.inl[0] + o = self.outl[0] + f = self.eta_volumetric_func + + if self.is_variable(i.p): + self.jacobian[k, i.p.J_col] = self.numeric_deriv(f, "p", i) + if self.is_variable(i.h): + self.jacobian[k, i.h.J_col] = self.numeric_deriv(f, "h", i) + if self.is_variable(o.p): + self.jacobian[k, o.p.J_col] = self.inl[0].calc_vol() + if self.is_variable(o.h): + self.jacobian[k, o.h.J_col] = -self.eta_vol.val + def flow_char_func(self): r""" Equation for given flow characteristic of a pump. diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index 2e0a06590..4cc1e1c3b 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -890,7 +890,7 @@ def fluid_balance_deriv(self, k, **kwargs): def calc_s(self): try: - return s_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=self.T.val_SI) + return s_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data, self.mixing_rule, T0=self.T.val_SI, solvent=self.solvent) except NotImplementedError: return np.nan @@ -911,18 +911,19 @@ def calc_results(self): if number_fluids > 1: h_from_T = h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule, solvent=self.solvent) if abs(h_from_T - self.h.val_SI) > ERR ** 0.5: - self.T.val_SI = np.nan - self.vol.val_SI = np.nan - self.v.val_SI = np.nan - self.s.val_SI = np.nan - msg = ( - "Could not find a feasible value for mixture temperature at " - f"connection {self.label}. The values for temperature, " - "specific volume, volumetric flow and entropy are set to nan. " - f"The deviation is {h_from_T - self.h.val_SI} J/kg." - ) - logger.error(msg) - _converged = False + # self.T.val_SI = np.nan + # self.vol.val_SI = np.nan + # self.v.val_SI = np.nan + # self.s.val_SI = np.nan + # msg = ( + # "Could not find a feasible value for mixture temperature at " + # f"connection {self.label}. The values for temperature, " + # "specific volume, volumetric flow and entropy are set to nan. " + # f"The deviation is {h_from_T - self.h.val_SI} J/kg." + # ) + # logger.error(msg) + # _converged = True + pass else: try: @@ -936,13 +937,11 @@ def calc_results(self): except ValueError: self.x.val_SI = np.nan - try: - if _converged: - self.vol.val_SI = self.calc_vol() - self.v.val_SI = self.vol.val_SI * self.m.val_SI - self.s.val_SI = self.calc_s() - except KeyError: - pass + if _converged: + self.vol.val_SI = self.calc_vol() + self.v.val_SI = self.vol.val_SI * self.m.val_SI + self.s.val_SI = self.calc_s() + for prop in fpd.keys(): self.get_attr(prop).val = convert_from_SI( prop, self.get_attr(prop).val_SI, self.get_attr(prop).unit diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index ede0271df..57772e481 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -2239,17 +2239,22 @@ def update_variables(self): relax = max(1, -2 * increment / container.val_SI) container.val_SI += increment / relax elif data["variable"] == "fluid": + if self.iter < 6: + # prevents fluid changes in first iteration! + relax = self.iter / 5 + else: + relax = 1 container = data["obj"].fluid container.val[data["fluid"]] += self.increment[ container.J_col[data["fluid"]] - ] + ] * relax if container.val[data["fluid"]] < ERR : container.val[data["fluid"]] = 0 elif container.val[data["fluid"]] > 1 - ERR : container.val[data["fluid"]] = 1 else: - # add increment + # component variables data["obj"].val += self.increment[data["obj"].J_col] # keep value within specified value range diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py index 8abb4bf35..1c10c6144 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/src/tespy/tools/fluid_properties/functions.py @@ -218,8 +218,8 @@ def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None, **kwargs): 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) + T = T_mix_ph(p, h , fluid_data, mixing_rule, T0, **kwargs) + return v_mix_pT(p, T, fluid_data, mixing_rule, **kwargs) def dv_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None): @@ -242,7 +242,7 @@ def v_mix_pT(p, T, fluid_data, mixing_rule=None, **kwargs): return 1 / pure_fluid["wrapper"].d_pT(p, T) else: _check_mixing_rule(mixing_rule, V_MIX_PT_DIRECT, "specific volume") - return V_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data) + return V_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data, **kwargs) def viscosity_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): diff --git a/src/tespy/tools/fluid_properties/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index f018fae0c..120d63535 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -241,6 +241,24 @@ def v_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return v +def v_mix_pT_incompressible_solution(p, T, fluid_data, **kwargs): + solvent = kwargs["solvent"] + + if p < fluid_data[solvent]["wrapper"].p_sat(T): + x_old = fluid_data[solvent]["mass_fraction"] + kwargs["x0"] = x_old + x_new = xsat_pT_incomp_solution(p, T, fluid_data, **kwargs) + x_water_gas = x_new - x_old + + v = 1 / fluid_data[solvent]["wrapper"].d_pT(p + 1e-6, T) * (1 - x_water_gas) + v += 1 / fluid_data["water"]["wrapper"].d_QT(1, T) * x_water_gas + fluid_data[solvent]["wrapper"].AS.set_mass_fractions([x_old]) + + return v + else: + return 1 / fluid_data[solvent]["wrapper"].d_pT(p, T) + + def viscosity_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): r""" Calculate dynamic viscosity from pressure and temperature. @@ -434,7 +452,8 @@ def cond_check(p, T, fluid_data, water_alias): 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, + "incomp-solution": v_mix_pT_incompressible_solution } diff --git a/src/tespy/tools/helpers.py b/src/tespy/tools/helpers.py index 019f306df..063e5711a 100644 --- a/src/tespy/tools/helpers.py +++ b/src/tespy/tools/helpers.py @@ -186,8 +186,7 @@ def latex_unit(unit): class UserDefinedEquation: - def __init__(self, label, func, deriv, conns, params={}, - latex={}): + def __init__(self, label, func, deriv, conns, params={}, latex={}): r""" A UserDefinedEquation allows use of generic user specified equations. @@ -610,10 +609,13 @@ def newton(func, deriv, params, y, **kwargs): valmin = kwargs.get('valmin', 70) valmax = kwargs.get('valmax', 3000) max_iter = kwargs.get('max_iter', 10) - tol_rel = kwargs.get('tol_rel', ERR ) - tol_abs = kwargs.get('tol_abs', ERR ) + tol_rel = kwargs.get('tol_rel', ERR) + tol_abs = kwargs.get('tol_abs', ERR ** 2) tol_mode = kwargs.get('tol_mode', 'abs') + if abs(y) <= 1e-6: + tol_mode = "abs" + # start newton loop expr = True i = 0 @@ -630,10 +632,12 @@ def newton(func, deriv, params, y, **kwargs): i += 1 if i > max_iter: - msg = ('Newton algorithm was not able to find a feasible value ' - 'for function ' + str(func) + '. Current value with x=' + - str(x) + ' is ' + str(func(params, x)) + - ', target value is ' + str(y) + '.') + msg = ( + 'The Newton algorithm was not able to find a feasible value ' + f'for function {func}. Current value with x={x} is ' + f'{func(params, x)}, target value is {y}, residual is {res} ' + f'after {i} iterations.' + ) logger.debug(msg) break @@ -649,7 +653,7 @@ def newton(func, deriv, params, y, **kwargs): def newton_with_kwargs( derivative, target_value, val0=300, valmin=70, valmax=3000, max_iter=10, - tol_rel=ERR, tol_abs=ERR, tol_mode="rel", **function_kwargs + tol_rel=ERR, tol_abs=ERR ** 2, tol_mode="rel", **function_kwargs ): # start newton loop @@ -660,6 +664,9 @@ def newton_with_kwargs( function = function_kwargs["function"] relax = 1 + if abs(target_value) <= 1e-6: + tol_mode = "abs" + while expr: # calculate function residual and new value function_kwargs[parameter] = x diff --git a/test_absorber.py b/test_absorber.py index 3425ec943..c876907f6 100644 --- a/test_absorber.py +++ b/test_absorber.py @@ -1,13 +1,9 @@ -from tespy.components.component import Component +from tespy.tools.data_containers import ComponentProperties as dc_cp +from sorption import Sorption -from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution - -class Absorber(Component): - - def __init__(self, label, **kwargs): - super().__init__(label, **kwargs) +class Absorber(Sorption): @staticmethod def inlets(): @@ -17,74 +13,15 @@ def inlets(): def outlets(): return ["out1"] - def preprocess(self, num_vars): - self.h2o = "water" - super().preprocess(num_vars) - def get_parameters(self): - return {} - - def get_mandatory_constraints(self): - constraints = { - 'mass_flow_constraints': { - 'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv, - 'constant_deriv': True,# 'latex': self.mass_flow_func_doc, - 'num_eq': 1}, - 'fluid_constraints': { - 'func': self.fluid_func, 'deriv': self.fluid_deriv, - 'constant_deriv': False,# 'latex': self.fluid_func_doc, - 'num_eq': 1}, - 'pressure_constraints': { - 'func': self.pressure_equality_func, - 'deriv': self.pressure_equality_deriv, - 'constant_deriv': True, - 'latex': self.pressure_equality_func_doc, - 'num_eq': 2}, - "saturation_constraints_libr": { - "func": self.saturated_solution_libr_func, - "deriv": self.saturated_solution_libr_deriv, - "constant_deriv": False, - "num_eq": 1 - }, + return { + 'Q': dc_cp( + max_val=0, + func=self.heat_func, + num_eq=1, + deriv=self.heat_deriv + ) } - if self.outl[0].solvent in self.outl[0].fluid.is_var: - constraints["saturation_constraints_water"] = { - "func": self.saturated_solution_water_func, - "deriv": self.saturated_solution_water_deriv, - "constant_deriv": False, - "num_eq": 1 - } - return constraints - - def mass_flow_func(self): - return self.inl[0].m.val_SI + self.inl[1].m.val_SI - self.outl[0].m.val_SI - - def mass_flow_deriv(self, k): - for c in self.inl: - if c.m.is_var: - self.jacobian[k, c.m.J_col] = 1 - - if self.outl[0].m.is_var: - self.jacobian[k, self.outl[0].m.J_col] = -1 - - def fluid_func(self): - return ( - self.inl[1].m.val_SI * self.inl[1].fluid.val["LiBr"] - - self.outl[0].m.val_SI * self.outl[0].fluid.val["LiBr"] - ) - - def fluid_deriv(self, increment_filter, k): - outl = self.outl[0] - inl = self.inl[1] - - if inl.m.is_var: - self.jacobian[k, inl.m.J_col] = inl.fluid.val[inl.solvent] - if inl.solvent in inl.fluid.is_var: - self.jacobian[k, inl.fluid.J_col[inl.solvent]] = inl.m.val_SI - if outl.m.is_var: - self.jacobian[k, outl.m.J_col] = -outl.fluid.val[outl.solvent] - if outl.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -outl.m.val_SI def pressure_equality_func(self): residual = [] @@ -108,37 +45,6 @@ def pressure_equality_deriv(self, k): self.jacobian[k, self.outl[0].p.J_col] = -1 k += 1 - def saturated_solution_water_func(self): - outl = self.outl[0] - return 1 - outl.fluid.val[outl.solvent] - outl.fluid.val[self.h2o] - - - def saturated_solution_water_deriv(self, increment_filter, k): - outl = self.outl[0] - if self.h2o in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[self.h2o]] = -1 - if outl.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -1 - - def saturated_solution_libr_func(self): - outl = self.outl[0] - x_previous = outl.fluid.val[outl.solvent] - T = outl.calc_T() - x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=outl.solvent, x0=x_previous) - outl.fluid_data[outl.solvent]["wrapper"].AS.set_mass_fractions([x_previous]) - return x_libr - outl.fluid.val[outl.solvent] - - def saturated_solution_libr_deriv(self, increment_filter, k): - outl = self.outl[0] - if outl.p.is_var: - deriv = self.numeric_deriv(self.saturated_solution_libr_func, "p", outl) - self.jacobian[k, outl.p.J_col] = deriv - if outl.h.is_var: - deriv = self.numeric_deriv(self.saturated_solution_libr_func, "h", outl) - self.jacobian[k, outl.h.J_col] = deriv - if outl.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[outl.solvent]] = self.numeric_deriv(self.saturated_solution_libr_func, outl.solvent, outl) - @staticmethod def is_branch_source(): return True @@ -154,11 +60,8 @@ def start_branch(self): return {outconn.label: branch} - def propagate_to_target(self, branch): - return - def propagate_wrapper_to_target(self, branch): - if branch["connections"][-1] == self.inl[0]: + if branch["connections"][-1] == self.inl[1]: branch["components"] += [self] return @@ -188,8 +91,8 @@ def propagate_wrapper_to_target(self, branch): absorber = Absorber("absorber") - c1 = Connection(water, "out1", absorber, "in1", label="1") - c2 = Connection(poor, "out1", absorber, "in2", label="2") + c1 = Connection(water, "out1", absorber, "in2", label="1") + c2 = Connection(poor, "out1", absorber, "in1", label="2") c3 = Connection(absorber, "out1", rich, "in1", label="3") nw.add_conns(c1, c2, c3) diff --git a/test_desorber.py b/test_desorber.py index 152f159d3..c1bfa3429 100644 --- a/test_desorber.py +++ b/test_desorber.py @@ -1,13 +1,8 @@ +from tespy.tools.data_containers import ComponentProperties as dc_cp +from sorption import Sorption -from tespy.components.component import Component -from tespy.tools.fluid_properties.mixtures import xsat_pT_incomp_solution - - -class Desorber(Component): - - def __init__(self, label, **kwargs): - super().__init__(label, **kwargs) +class Desorber(Sorption): @staticmethod def inlets(): @@ -17,77 +12,15 @@ def inlets(): def outlets(): return ["out1", "out2"] - def preprocess(self, num_vars): - self.h2o = "water" - super().preprocess(num_vars) - def get_parameters(self): - return {} - - def get_mandatory_constraints(self): - constraints = { - 'mass_flow_constraints': { - 'func': self.mass_flow_func, 'deriv': self.mass_flow_deriv, - 'constant_deriv': True,# 'latex': self.mass_flow_func_doc, - 'num_eq': 1}, - 'fluid_constraints': { - 'func': self.fluid_func, 'deriv': self.fluid_deriv, - 'constant_deriv': False,# 'latex': self.fluid_func_doc, - 'num_eq': 1}, - 'pressure_constraints': { - 'func': self.pressure_equality_func, - 'deriv': self.pressure_equality_deriv, - 'constant_deriv': True, - 'latex': self.pressure_equality_func_doc, - 'num_eq': 2}, - "saturation_constraints_libr": { - "func": self.saturated_solution_libr_func, - "deriv": self.saturated_solution_libr_deriv, - "constant_deriv": False, - "num_eq": 1 - }, + return { + 'Q': dc_cp( + min_val=0, + func=self.heat_func, + num_eq=1, + deriv=self.heat_deriv + ) } - if self.outl[1].solvent in self.outl[1].fluid.is_var: - constraints["saturation_constraints_water"] = { - "func": self.saturated_solution_water_func, - "deriv": self.saturated_solution_water_deriv, - "constant_deriv": False, - "num_eq": 1 - } - return constraints - - def mass_flow_func(self): - return self.inl[0].m.val_SI - self.outl[0].m.val_SI - self.outl[1].m.val_SI - - def mass_flow_deriv(self, k): - inl = self.inl[0] - if inl.m.is_var: - self.jacobian[k, inl.m.J_col] = 1 - - for outl in self.outl: - if outl.m.is_var: - self.jacobian[k, outl.m.J_col] = -1 - - def fluid_func(self): - inl = self.inl[0] - outl = self.outl[1] - return ( - inl.m.val_SI * inl.fluid.val[inl.solvent] - - outl.m.val_SI * outl.fluid.val[outl.solvent] - ) - - def fluid_deriv(self, increment_filter, k): - outl = self.outl[1] - inl = self.inl[0] - - if inl.m.is_var: - self.jacobian[k, inl.m.J_col] = inl.fluid.val[inl.solvent] - if inl.solvent in inl.fluid.is_var: - self.jacobian[k, inl.fluid.J_col[inl.solvent]] = inl.m.val_SI - if outl.m.is_var: - self.jacobian[k, outl.m.J_col] = -outl.fluid.val[outl.solvent] - if outl.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -outl.m.val_SI def pressure_equality_func(self): residual = [] @@ -111,36 +44,6 @@ def pressure_equality_deriv(self, k): self.jacobian[k, self.inl[0].p.J_col] = -1 k += 1 - def saturated_solution_water_func(self): - outl = self.outl[1] - return 1 - outl.fluid.val[outl.solvent] - outl.fluid.val[self.h2o] - - def saturated_solution_water_deriv(self, increment_filter, k): - outl = self.outl[1] - if self.h2o in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[self.h2o]] = -1 - if outl.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[outl.solvent]] = -1 - - def saturated_solution_libr_func(self): - outl = self.outl[1] - x_previous = outl.fluid.val[outl.solvent] - T = outl.calc_T() - x_libr = xsat_pT_incomp_solution(outl.p.val_SI, T, outl.fluid_data, solvent=outl.solvent, x0=x_previous) - outl.fluid_data[outl.solvent]["wrapper"].AS.set_mass_fractions([x_previous]) - return x_libr - outl.fluid.val[outl.solvent] - - def saturated_solution_libr_deriv(self, increment_filter, k): - outl = self.outl[1] - if outl.p.is_var: - deriv = self.numeric_deriv(self.saturated_solution_libr_func, "p", outl) - self.jacobian[k, outl.p.J_col] = deriv - if outl.h.is_var: - deriv = self.numeric_deriv(self.saturated_solution_libr_func, "h", outl) - self.jacobian[k, outl.h.J_col] = deriv - if outl.solvent in outl.fluid.is_var: - self.jacobian[k, outl.fluid.J_col[outl.solvent]] = self.numeric_deriv(self.saturated_solution_libr_func, outl.solvent, outl) - @staticmethod def is_branch_source(): return True @@ -174,14 +77,11 @@ def start_fluid_wrapper_branch(self): return branches - def propagate_to_target(self, branch): - return - def propagate_wrapper_to_target(self, branch): if self in branch["components"]: return - outconn = self.outl[1] + outconn = self.outl[0] branch["connections"] += [outconn] branch["components"] += [self] outconn.target.propagate_wrapper_to_target(branch) @@ -207,8 +107,8 @@ def propagate_wrapper_to_target(self, branch): valve = SimpleHeatExchanger("valve") c5 = Connection(rich, "out1", desorber, "in1", label="5") - c10 = Connection(desorber, "out1", water, "in1", label="10") - c6 = Connection(desorber, "out2", valve, "in1", label="6") + c10 = Connection(desorber, "out2", water, "in1", label="10") + c6 = Connection(desorber, "out1", valve, "in1", label="6") c7 = Connection(valve, "out1", poor, "in1", label="7") nw.add_conns(c5, c10, c6, c7) From caf52f7699ff3cba626b0cb2cc70b867a0b6e65f Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Wed, 14 Feb 2024 08:37:47 +0100 Subject: [PATCH 21/22] Move files to individual repository --- sorption_reference.py | 116 -------------- src/tespy/networks/network.py | 5 +- src/tespy/tools/fluid_properties/wrappers.py | 9 -- test_absorber.py | 154 ------------------ test_desorber.py | 157 ------------------- 5 files changed, 2 insertions(+), 439 deletions(-) delete mode 100644 sorption_reference.py delete mode 100644 test_absorber.py delete mode 100644 test_desorber.py diff --git a/sorption_reference.py b/sorption_reference.py deleted file mode 100644 index 616e084d0..000000000 --- a/sorption_reference.py +++ /dev/null @@ -1,116 +0,0 @@ - -from CoolProp.CoolProp import AbstractState -import CoolProp as CP -import numpy as np - - -libr = AbstractState("INCOMP", "LiBr") - -libr.set_mass_fractions([0.5]) -libr.update(CP.QT_INPUTS, 0, 300) -libr.p() -libr.set_mass_fractions([0.2]) -libr.update(CP.QT_INPUTS, 0, 300) -libr.p() - -# saturation pressure by temperature and libr mass fraction x - -def psat_TX(T, X): - libr.set_mass_fractions([X]) - libr.update(CP.QT_INPUTS, 0, T) - return libr.p() - - -def Xsat_Tp(T, p): - X_min = 0 - X_max = 0.75 - X = .3 - d = 1e-5 - while True: - res = psat_TX(T, X) - p - deriv = (psat_TX(T, X + d) - psat_TX(T, X - d)) / (2 * d) - X -= res / deriv - # print(X) - - if abs(res) < 1e-6: - if X < X_min: - return X_min - elif X > X_max: - return X_max - break - - if X >= X_max: - X = X_max - d - elif X <= X_min: - X = X_min + d - - psat_TX(T, X) - return X - - -def Tsat_pX(p, X): - return newton(psat_TX, p, X) - - -def newton(func, p, X): - d = 1e-3 - # T = 300 - T_min = libr.trivial_keyed_output(CP.iT_min) - T_max = libr.trivial_keyed_output(CP.iT_max) - T = (T_min + T_max) / 2 - while True: - res = func(T, X) - p - deriv = (func(T + d, X) - func(T -d, X)) / (2 * d) - T -= res / deriv - - if T >= T_max: - T = T_max - d - elif T <= T_min: - T = T_min + d - - if abs(res) < 1e-6: - return T - -# watercycle -T_evap = 275.15 -p_evap = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_evap, "water") - -T_cond = 309.15 -p_cond = CP.CoolProp.PropsSI("P", "Q", 0, "T", T_cond, "water") - -# heat source temperature -T_sol_des_out = 350 - -# heat sink temperature -T_sol_abs_out = 306.15 - -x_rich = Xsat_Tp(T_sol_abs_out, p_evap) - -x_water_rich = 1 - x_rich -T_water_abs_in = T_evap -h_water_abs_in = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_abs_in, "water") - -h_sol_abs_out = libr.hmass() -s_sol_abs_out = libr.smass() - -eta_s_pump = 0.9 -libr.update(CP.PSmass_INPUTS, p_cond, s_sol_abs_out) -h_pump_out_s = libr.hmass() -h_pump_out = h_sol_abs_out + (h_pump_out_s - h_sol_abs_out) / eta_s_pump - -x_poor = Xsat_Tp(T_sol_des_out, p_cond) -h_poor = libr.hmass()#update(CP.PT_INPUTS, p_cond, T_sol_des_out) - -delta_T = 0 -T_water_des_out = T_sol_des_out - delta_T -h_water_des_out = CP.CoolProp.PropsSI("H", "Q", 1, "T", T_water_des_out, "water") - -x_water_poor = 1 - x_poor - -m_water = 1 -m_rich = m_water * (1 - x_water_poor) / (x_water_rich - x_water_poor) - -m_poor = m_rich - m_water - -delta_H_abs = -m_water * h_water_abs_in - m_poor * h_poor + m_rich * h_sol_abs_out -delta_H_des = m_water * h_water_des_out + m_poor * h_poor - m_rich * h_pump_out diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index 57772e481..95803bb86 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -925,8 +925,8 @@ def propagate_fluid_wrappers(self): num_potential_fluids = len(potential_fluids) if num_potential_fluids == 0: msg = ( - "The follwing connections of your network are missing any " - "kind of fluid composition information:" + "The following connections of your network are missing any " + "kind of fluid composition information: " + ", ".join([c.label for c in all_connections]) + "." ) raise hlp.TESPyNetworkError(msg) @@ -2228,7 +2228,6 @@ def matrix_inversion(self): self.increment = self.residual * 0 def update_variables(self): - # add the increment for data in self.variables_dict.values(): if data["variable"] in ["m", "h"]: container = data["obj"].get_attr(data["variable"]) diff --git a/src/tespy/tools/fluid_properties/wrappers.py b/src/tespy/tools/fluid_properties/wrappers.py index 01099033e..90c5b2034 100644 --- a/src/tespy/tools/fluid_properties/wrappers.py +++ b/src/tespy/tools/fluid_properties/wrappers.py @@ -186,15 +186,6 @@ def h_ps(self, p, s): return self.AS.hmass() def h_pT(self, p, T): - if self.back_end == "INCOMP": - if T == (self._T_max + self._T_min) / 2: - T += ERR - self.AS.update(CP.PT_INPUTS, p, T) - h = self.AS.hmass() * 0.5 - T -= 2 * ERR - self.AS.update(CP.PT_INPUTS, p, T) - h += self.AS.hmass() * 0.5 - return h self.AS.update(CP.PT_INPUTS, p, T) return self.AS.hmass() diff --git a/test_absorber.py b/test_absorber.py deleted file mode 100644 index c876907f6..000000000 --- a/test_absorber.py +++ /dev/null @@ -1,154 +0,0 @@ - -from tespy.tools.data_containers import ComponentProperties as dc_cp -from sorption import Sorption - - -class Absorber(Sorption): - - @staticmethod - def inlets(): - return ["in1", "in2"] - - @staticmethod - def outlets(): - return ["out1"] - - def get_parameters(self): - return { - 'Q': dc_cp( - max_val=0, - func=self.heat_func, - num_eq=1, - deriv=self.heat_deriv - ) - } - - def pressure_equality_func(self): - residual = [] - for c in self.inl : - residual += [c.p.val_SI - self.outl[0].p.val_SI] - return residual - - def pressure_equality_deriv(self, k): - r""" - Calculate partial derivatives for all pressure equations. - - Returns - ------- - deriv : ndarray - Matrix with partial derivatives for the fluid equations. - """ - for c in self.inl: - if c.p.is_var: - self.jacobian[k, c.p.J_col] = 1 - if self.outl[0].p.is_var: - self.jacobian[k, self.outl[0].p.J_col] = -1 - k += 1 - - @staticmethod - def is_branch_source(): - return True - - def start_branch(self): - outconn = self.outl[0] - branch = { - "connections": [outconn], - "components": [self, outconn.target], - "subbranches": {} - } - outconn.target.propagate_to_target(branch) - - return {outconn.label: branch} - - def propagate_wrapper_to_target(self, branch): - if branch["connections"][-1] == self.inl[1]: - branch["components"] += [self] - return - - if self in branch["components"]: - return - - outconn = self.outl[0] - branch["connections"] += [outconn] - branch["components"] += [self] - outconn.target.propagate_wrapper_to_target(branch) - - -if __name__ == "__main__": - - from tespy.components import Source, Sink - from tespy.networks import Network - from tespy.connections import Connection - - # this runs the references calculations and imports all variables from there - from sorption_reference import * - - nw = Network() - - water = Source("water source") - rich = Sink("rich solution") - poor = Source("poor solution") - - absorber = Absorber("absorber") - - c1 = Connection(water, "out1", absorber, "in2", label="1") - c2 = Connection(poor, "out1", absorber, "in1", label="2") - c3 = Connection(absorber, "out1", rich, "in1", label="3") - - nw.add_conns(c1, c2, c3) - - c1.set_attr(fluid={"water": 1}, p=0.01e5, m=1, x=1) - c2.set_attr(fluid={"water": x_water_poor, "INCOMP::LiBr": 1 - x_water_poor}, h=h_poor, mixing_rule="incomp-solution", solvent="LiBr") - c3.set_attr(h=h_sol_abs_out, p0=0.01e5, m0=m_rich) - - nw.solve("design") - nw.print_results() - print(c3.fluid.val) - - c2.set_attr(m=10) - c3.set_attr(h=None) - - nw.solve("design") - nw.print_results() - - print(c1.m.val_SI * c1.h.val_SI + c2.m.val_SI * c2.h.val_SI - c3.m.val_SI * c3.h.val_SI) - - # how does the reference temperature affect the energy balance - # eb_prev = 0 - # p_ref = 1e5 - # for T_ref in range(274, 350): - # h_water_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, "water") - # h_poor_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{c2.fluid.val['LiBr']}]") - # h_rich_ref = CP.CoolProp.PropsSI("H", "P", p_ref, "T", T_ref, f"INCOMP::LiBr[{c3.fluid.val['LiBr']}]") - # (CP.CoolProp.PropsSI("Cpmass", "P", p_ref, "T", T_ref, "water")) - # eb = (c1.h.val_SI - h_water_ref) * c1.m.val_SI + (c2.h.val_SI - h_poor_ref) * c2.m.val_SI - (c3.h.val_SI - h_rich_ref) * c3.m.val_SI - # eb_prev = eb - - # s - - # some checks - # nw.print_results() - print(m_rich, c3.m.val_SI) - assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) - assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) - - # use temperature instead of enthalpy - c2.set_attr(m=None) - c3.set_attr(T=306.15) - nw.solve("design") - nw.print_results() - # check results, should be identical to previous ones - assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) - assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) - - # fix the either the water or the libr mass fraction and the temperature or pressure - # gives us the other respective value, pressure in this case - c1.set_attr(p=None) - c3.set_attr(fluid={"INCOMP::LiBr": 0.50}, T=305) - nw.solve("design") - assert round(c1.m.val_SI * c1.fluid.val["water"] + c2.m.val_SI * c2.fluid.val["water"], 4) == round(c3.m.val_SI * c3.fluid.val["water"], 4) - assert round(c2.m.val_SI * c2.fluid.val["LiBr"], 4) == round(c3.m.val_SI * c3.fluid.val["LiBr"], 4) - - print(f"T3: {c3.T.val_SI}, p3: {c3.p.val_SI}, fluid3: {c3.fluid.val}") - - nw.print_results() diff --git a/test_desorber.py b/test_desorber.py deleted file mode 100644 index c1bfa3429..000000000 --- a/test_desorber.py +++ /dev/null @@ -1,157 +0,0 @@ -from tespy.tools.data_containers import ComponentProperties as dc_cp -from sorption import Sorption - - -class Desorber(Sorption): - - @staticmethod - def inlets(): - return ["in1"] - - @staticmethod - def outlets(): - return ["out1", "out2"] - - def get_parameters(self): - return { - 'Q': dc_cp( - min_val=0, - func=self.heat_func, - num_eq=1, - deriv=self.heat_deriv - ) - } - - def pressure_equality_func(self): - residual = [] - for c in self.outl: - residual += [c.p.val_SI - self.inl[0].p.val_SI] - return residual - - def pressure_equality_deriv(self, k): - r""" - Calculate partial derivatives for all pressure equations. - - Returns - ------- - deriv : ndarray - Matrix with partial derivatives for the fluid equations. - """ - for c in self.outl: - if c.p.is_var: - self.jacobian[k, c.p.J_col] = 1 - if self.inl[0].p.is_var: - self.jacobian[k, self.inl[0].p.J_col] = -1 - k += 1 - - @staticmethod - def is_branch_source(): - return True - - @staticmethod - def is_wrapper_branch_source(): - return True - - def start_branch(self): - branches = {} - for outconn in self.outl: - branch = { - "connections": [outconn], - "components": [self, outconn.target], - "subbranches": {} - } - outconn.target.propagate_to_target(branch) - - branches[outconn.label] = branch - return branches - - def start_fluid_wrapper_branch(self): - branches = {} - for outconn in self.outl: - branch = { - "connections": [outconn], - "components": [self] - } - outconn.target.propagate_wrapper_to_target(branch) - branches[outconn.label] = branch - - return branches - - def propagate_wrapper_to_target(self, branch): - if self in branch["components"]: - return - - outconn = self.outl[0] - branch["connections"] += [outconn] - branch["components"] += [self] - outconn.target.propagate_wrapper_to_target(branch) - - -if __name__ == "__main__": - - from tespy.components import Source, Sink, SimpleHeatExchanger - from tespy.networks import Network - from tespy.connections import Connection - - # this runs the references calculations and imports all variables from there - from sorption_reference import * - - - nw = Network() - - water = Sink("water source") - rich = Source("rich solution") - poor = Sink("poor solution") - - desorber = Desorber("desorber") - valve = SimpleHeatExchanger("valve") - - c5 = Connection(rich, "out1", desorber, "in1", label="5") - c10 = Connection(desorber, "out2", water, "in1", label="10") - c6 = Connection(desorber, "out1", valve, "in1", label="6") - c7 = Connection(valve, "out1", poor, "in1", label="7") - - nw.add_conns(c5, c10, c6, c7) - - c5.set_attr(fluid={"water": x_water_rich, "INCOMP::LiBr": 1 - x_water_rich}, p=p_cond, m=10, h=h_pump_out, mixing_rule="incomp-solution", solvent="LiBr") - c10.set_attr(fluid={"water": 1}, x=1) - c6.set_attr(h=h_poor, p0=0.1e5, m0=m_rich) - c7.set_attr(p=p_cond, h=h_poor * 1.2) - - nw.solve("design") - print(c6.fluid.val) - print(c10.h.val_SI * c10.m.val_SI + c6.h.val_SI * c6.m.val_SI - c5.m.val_SI * c5.h.val_SI) - print(c7.h.val_SI) - print(c7.T.val_SI) - - c10.set_attr(x=None, T=320) - nw.solve("design") - print(c6.fluid.val) - print(c10.h.val_SI * c10.m.val_SI + c6.h.val_SI * c6.m.val_SI - c5.m.val_SI * c5.h.val_SI) - exit() - nw.print_results() - - print(T_sol_des_out) - print(h_poor) - print(c6.h.val_SI) - print(c6.T.val_SI) - print(c10.T.val_SI) - print(c6.fluid.val) - c6.set_attr(m=9) - c6.set_attr(h=None) - - nw.solve("design") - nw.print_results() - print(c6.fluid.val) - - c6.set_attr(m=None, T=350) - nw.solve("design") - nw.print_results() - print(c6.fluid.val) - - c6.set_attr(fluid={"INCOMP::LiBr": 0.7}) - c5.set_attr(p=None) - - nw.solve("design") - nw.print_results() - print(c6.fluid.val) From 6e48930c47c392b556724a33b0948caa4e48bc00 Mon Sep 17 00:00:00 2001 From: Francesco Witte Date: Thu, 15 Feb 2024 17:21:05 +0100 Subject: [PATCH 22/22] Allow specific exergy calculations for incompressible solutions --- src/tespy/connections/connection.py | 3 ++- src/tespy/tools/fluid_properties/functions.py | 10 +++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index 4cc1e1c3b..1edb8aaf9 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -1086,7 +1086,8 @@ def get_physical_exergy(self, pamb, Tamb): """ self.ex_therm, self.ex_mech = fp.functions.calc_physical_exergy( self.h.val_SI, self.s.val_SI, self.p.val_SI, pamb, Tamb, - self.fluid_data, self.mixing_rule, self.T.val_SI + self.fluid_data, self.mixing_rule, self.T.val_SI, + solvent=self.solvent ) self.Ex_therm = self.ex_therm * self.m.val_SI self.Ex_mech = self.ex_mech * self.m.val_SI diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py index 1c10c6144..c9a0a5ba3 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/src/tespy/tools/fluid_properties/functions.py @@ -34,7 +34,7 @@ def isentropic(p_1, h_1, p_2, fluid_data, mixing_rule=None, T0=None, **kwargs): return h_mix_pT(p_2, T_2, fluid_data, mixing_rule, **kwargs) -def calc_physical_exergy(h, s, p, pamb, Tamb, fluid_data, mixing_rule=None, T0=None): +def calc_physical_exergy(h, s, p, pamb, Tamb, fluid_data, mixing_rule=None, T0=None, **kwargs): r""" Calculate specific physical exergy. @@ -65,11 +65,11 @@ def calc_physical_exergy(h, s, p, pamb, Tamb, fluid_data, mixing_rule=None, T0=N e^\mathrm{PH} = e^\mathrm{T} + e^\mathrm{M} """ - h_T0_p = h_mix_pT(p, Tamb, fluid_data, mixing_rule) - s_T0_p = s_mix_pT(p, Tamb, fluid_data, mixing_rule) + h_T0_p = h_mix_pT(p, Tamb, fluid_data, mixing_rule, **kwargs) + s_T0_p = s_mix_pT(p, Tamb, fluid_data, mixing_rule, **kwargs) ex_therm = (h - h_T0_p) - Tamb * (s - s_T0_p) - h0 = h_mix_pT(pamb, Tamb, fluid_data, mixing_rule) - s0 = s_mix_pT(pamb, Tamb, fluid_data, mixing_rule) + h0 = h_mix_pT(pamb, Tamb, fluid_data, mixing_rule, **kwargs) + s0 = s_mix_pT(pamb, Tamb, fluid_data, mixing_rule, **kwargs) ex_mech = (h_T0_p - h0) - Tamb * (s_T0_p - s0) return ex_therm, ex_mech