Add sorption features to TESPy

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 = {

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 = {

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)

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:

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:]:

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,
@@ -211,7 +216,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
             )
         )
@@ -344,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.
@@ -511,7 +537,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)
 

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:
@@ -719,6 +720,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[self.solvent]["wrapper"].AS.set_mass_fractions([self.fluid.val[self.solvent]])
 
     def primary_ref_func(self, k, **kwargs):
         variable = kwargs["variable"]
@@ -739,23 +742,23 @@ 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
 
     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):
@@ -770,28 +773,28 @@ 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):
         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
 
@@ -887,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
 
@@ -902,23 +905,25 @@ def solve(self, increment_filter):
             data.deriv(k, **data.func_params)
 
     def calc_results(self):
+        _converged = True
         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
+            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 = True
+                pass
 
         else:
             try:
@@ -1023,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
@@ -1081,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

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:
@@ -893,16 +891,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."
@@ -913,14 +925,15 @@ 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)
 
             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]
@@ -1761,7 +1774,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
 
@@ -2215,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"])
@@ -2226,17 +2238,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