Implement an IncompressibleFluidWrapper class

docs/advanced_features/fluid_properties.rst

@@ -4,8 +4,13 @@ Fluid properties
 ================
 The default fluid property engine `CoolProp <https://coolprop.org/>`_. All
 available fluids can be found on their homepage. Also see :cite:`Bell2014`.
-Since version 0.7 of TESPy it is possible to use other engines. TESPy comes with
-two additional predefined engines, i.e.
+Since version 0.7 of TESPy it is possible to use other engines. TESPy supports
+one additional predefined engine, which can be used to fit functions to simple
+custom fluid property data, the :code:`IncompressibleFluidWrapper`. See
+:ref:`this section <incompressible_wrapper_label>` for more information.
+
+On top, there are two additional predefined engines, which are untested but may
+serve as an inspiration for you to create your onw one, i.e.
 
 - the `iapws <https://github.com/jjgomera/iapws/>`_ library and
 - the `pyromat <https://github.com/chmarti1/PYroMat/>`_ library.
@@ -62,19 +67,131 @@ If you are looking for heat transfer fluids, the list of incompressible
 might be interesting for you. In contrast to the pure fluids, the properties
 cover liquid state only.
 
+If you have measurement data or data from manufacturer data sheets and want to
+use these in your TESPy model, TESPy can create fitting functions for these
+with the :code:`IncompressibleFluidWrapper`. See
+:ref:`this section <incompressible_wrapper_label>` for more information.
+
 Fluid mixtures
 ++++++++++++++
 TESPy provides support for three types of mixtures:
 
 - ideal: Mixtures for gases only.
 - ideal-cond: Mixture for gases with condensation calculation for water share.
-- incompressible: Mixtures for incompressible fluids.
+- incompressible: Mixtures for CoolProp-based incompressible fluids.
+
+These mixtures are handled externally by TESPy by using the pure fluid
+properties of CoolProp and then applying the respective mixing rules, read more
+about it :ref:`here <mixture_routines_label>`.
+
+More accurate formulations are available directly through CoolProp, which
+provides a back end for predefined mixtures. This back end is rather instable
+when using HEOS. In general, to use the mixture feature of CoolProp we
+recommend using the REFPROP back end instead of HEOS. Also note, the CoolProp
+mixture back end is not tested thoroughly. Please reach out if you would like
+to support us in adopting the TESPy implementation.
+
+.. _incompressible_wrapper_label:
+
+IncompressibleFluidWrapper
+--------------------------
+You can use the :code:`IncompressibleFluidWrapper` engine of TESPy to model a
+fluid based on your own data. To do this, you need tabular data as function of
+temperature:
+
+- Mass density
+- Mass specific heat capacity
+- Dynamic viscosity
+
+The :code:`IncompressibleFluidWrapper` will automatically fit functions to your
+data:
+
+- density and heat capacity through linear interpolation:
+  :math:`f\left(T\right) = A + B \cdot T`
+- viscosity through an exponential polynomial equation
+  :math:`\eta\left(T\right) = e ^ {\frac{A}{T ^ 3} + \frac{B}{T ^ 2} + \frac{C}{T} + D}`
+
+We can make use of the engine as shown in the example below. First, we set up
+a very simple system, just a flow of fluid through a heat exchanger.
+
+.. code-block:: python
+
+    >>> from tespy.components import Sink
+    >>> from tespy.components import Source
+    >>> from tespy.components import SimpleHeatExchanger
+    >>> from tespy.connections import Connection
+    >>> from tespy.networks import Network
+    >>> from tespy.tools.fluid_properties import IncompressibleFluidWrapper
+    >>> import numpy as np
+    >>> nw = Network(iterinfo=False)
+    >>> nw.units.set_defaults(
+    ...     temperature="°C",
+    ...     pressure="bar",
+    ...     heat="kW"
+    ... )
+
+    >>> heatexchanger = SimpleHeatExchanger("heat exchanger")
+
+    >>> so = Source("source")
+    >>> si = Sink("sink")
+
+    >>> c1 = Connection(so, "out1", heatexchanger, "in1", label="c1")
+    >>> c2 = Connection(heatexchanger, "out1", si, "in1", label="c2")
+
+    >>> nw.add_conns(c1, c2)
 
-Furthermore, CoolProp provides a back end for predefined mixtures, which is
-rather instable using HEOS. Using the CoolProp mixture back-end is not tested,
-reach out if you would like to support us in adopting the TESPy implementation.
-In general, to use the mixture feature of CoolProp we recommend using the
-REFPROP back end instead of HEOS.
+Next, we have to prepare our data to be utilized by TESPy. The information
+required needs to be passed in SI units and must be gridded with the same
+spacing of temperature for all measurements.
+
+.. attention
+
+    Please note, that this example is purely for showing how to utilize this
+    implementation. In the example only 2 datapoints are available. This works
+    well for density and heat capacity as linear fitting functions are applied.
+    For viscosity this is somewhat sketchy, since you would need at least 4
+    datapoints to fit a function of polynomial 3!
+
+.. code-block:: python
+
+    >>> fluid_kwargs = {
+    ...     "temperature_data": np.array([273.15, 373.15]),  # K
+    ...     "density_data": np.array([1000, 1100]),  # kg/m3
+    ...     "heat_capacity_data": np.array([4000, 4100]) * 1e3,  # J/kg
+    ...     "viscosity_data": np.array([0.05, 0.00025])  # Pa*s
+    ... }
+
+.. attention::
+
+    The keys of the dictionary must be named :code:`temperature_data`,
+    :code:`density_data`, :code:`heat_capacity_data` and
+    :code:`viscosity_data`!
+
+Then, we can specify the fluid on one of our connections by providing two
+additional keywords:
+
+- :code:`fluid_engines` indicating which fluid uses which wrapper class.
+- :code:`fluid_wrapper_kwargs` providing for each fluid (if applicable) the
+  required data.
+
+.. code-block:: python
+
+    >>> c1.set_attr(
+    ...     fluid={"f": 1},
+    ...     fluid_engines={"f": IncompressibleFluidWrapper},
+    ...     fluid_wrapper_kwargs={"f": fluid_kwargs},
+    ... )
+
+We can specify missing boundary conditions and solve the problem.
+
+.. code-block:: python
+
+    >>> c1.set_attr(v=2, p=1, T=30)
+    >>> c2.set_attr(p=0.9)
+    >>> heatexchanger.set_attr(Q=10)
+
+    >>> nw.solve("design")
+    >>> nw.assert_convergence()
 
 Using other engines
 -------------------
@@ -150,6 +267,11 @@ enthalpy to the temperature. Lastly, to make the calculation of isentropic
 efficiencies possible, we can add the equation for change in enthalpy on
 isentropic change of pressure for an ideal gas.
 
+.. tip::
+
+    You can also inject :code:`kwargs` from the connection specification to
+    your concrete wrapper instance, see the section on
+    :ref:`incompressible fluids <incompressible_wrapper_label>`.
 
 .. code-block:: python
 
@@ -160,7 +282,7 @@ isentropic change of pressure for an ideal gas.
 
     >>> class KKHWrapper(FluidPropertyWrapper):
     ...
-    ...     def __init__(self, fluid, back_end=None, reference_temperature=298.15) -> None:
+    ...     def __init__(self, fluid, back_end=None, reference_temperature=298.15, **kwargs) -> None:
     ...         super().__init__(fluid, back_end)
     ...
     ...         if self.fluid not in COEF:
@@ -271,7 +393,6 @@ the previous section:
 
 .. code-block:: python
 
-
     >>> from tespy.components import Sink
     >>> from tespy.components import Source
     >>> from tespy.components import Turbine
@@ -304,6 +425,8 @@ the previous section:
     >>> round(c2.T.val, 1)
     306.3
 
+.. _mixture_routines_label:
+
 Mixture routines in TESPy
 -------------------------
 Different types of mixture routines are implemented in TESPy. You can select,

docs/knowledge_center/faq.rst

@@ -233,6 +233,19 @@ Customizing model behavior
             ...         for var in [c.p, c.h]
             ...     ]
 
+    .. dropdown:: I have custom fluid property data, how can I integrate them into my model?
+
+        In the documentation section
+        :ref:`on fluid property engines <fluid_properties_label>` you will find
+        a lot of helpful information. Specifically for liquids/incompressibles
+        there is already an interface, that takes your datatables and then fits
+        equations to them and integrates them  into your tespy model. Check out
+        the information on that topic in
+        :ref:`this section <incompressible_wrapper_label>`. You can also
+        implement your own class, that handles your fluid property equations.
+        Follow the examples given in the sections mentioned to learn, how that
+        can be accomplished.
+
 .. _faq_postprocessing_label:
 
 Visualization, post-processing and cycle analysis

docs/whats_new/v0-9-12.rst

@@ -1,6 +1,28 @@
 v0.9.12 - Under development
 +++++++++++++++++++++++++++
 
+New Features
+############
+- Custom fluid property wrappers can now receive arbitrary :code:`kwargs`
+  making injection of information for the underlying property models much
+  easier. Check the section on the
+  :ref:`incompressible fluid properties <incompressible_wrapper_label>`
+ (`PR #877 <https://github.com/oemof/tespy/pull/877>`__).
+- There is a new :code:`FluidPropertyWrapper` for incompressible fluids such as
+  thermo-oils, which you can pass measurement or manufacturer data to. The
+  :code:`IncompressibleFluidWrapper` will automatically fit functions to the
+  data points you provide:
+
+  - heat capacity and density: linear interpolation
+    :math:`f\left(T\right) = A + B \cdot T`.
+  - viscosity: exponential polynomial equation
+    :math:`\eta\left(T\right) = e ^ {\frac{A}{T ^ 3} + \frac{B}{T ^ 2} + \frac{C}{T} + D}`
+
+  For an example in a TESPy model see
+  :ref:`this section <incompressible_wrapper_label>`
+
+  (`PR #878 <https://github.com/oemof/tespy/pull/878>`__)
+
 Other changes
 #############
 - The optimization API has changed to integrate :code:`pymoo` instead of
@@ -63,7 +85,6 @@ Bug Fixes
   pressure relaxation factor to overwrite the factors for all other variables
   (`PR #875 <https://github.com/oemof/tespy/pull/875>`__).
 
-
 Contributors
 ############
 - Francesco Witte (`@fwitte <https://github.com/fwitte>`__)

src/tespy/connections/connection.py

@@ -772,6 +772,9 @@ def _fluid_specification(self, key, value):
         elif key == "fluid_balance":
             self.fluid_balance.is_set = value
 
+        elif key == "fluid_wrapper_kwargs":
+            self.fluid.wrapper_kwargs = value
+
         else:
             msg = f"Connections do not have an attribute named {key}"
             logger.error(msg)
@@ -833,6 +836,7 @@ def _create_fluid_wrapper(self):
         for fluid in self.fluid.val:
             if fluid in self.fluid.wrapper:
                 continue
+
             if fluid not in self.fluid.engine:
                 self.fluid.engine[fluid] = CoolPropWrapper
 
@@ -842,7 +846,13 @@ def _create_fluid_wrapper(self):
             else:
                 self.fluid.back_end[fluid] = None
 
-            self.fluid.wrapper[fluid] = self.fluid.engine[fluid](fluid, back_end)
+            wrapper_kwargs = {}
+            if fluid in self.fluid.wrapper_kwargs:
+                wrapper_kwargs = self.fluid.wrapper_kwargs[fluid]
+
+            self.fluid.wrapper[fluid] = self.fluid.engine[fluid](
+                fluid, back_end, **wrapper_kwargs
+            )
 
     def _precalc_guess_values(self):
         """

src/tespy/networks/network.py

@@ -917,16 +917,21 @@ def _propagate_fluid_wrappers(self):
             any_fluids_set = []
             engines = {}
             back_ends = {}
+            wrapper_kwargs = {}
             any_fluids = []
             any_fluids0 = []
             mixing_rules = []
             for c in all_connections:
                 for f in c.fluid.is_set:
                     any_fluids_set += [f]
+
                     if f in c.fluid.engine:
                         engines[f] = c.fluid.engine[f]
                     if f in c.fluid.back_end:
                         back_ends[f] = c.fluid.back_end[f]
+                    if f in c.fluid.wrapper_kwargs:
+                        wrapper_kwargs[f] = c.fluid.wrapper_kwargs[f]
+
                 any_fluids += list(c.fluid.val.keys())
                 any_fluids0 += list(c.fluid.val0.keys())
                 if c.mixing_rule is not None:
@@ -957,7 +962,7 @@ 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 "
+                    "The following connections of your network are missing any "
                     "kind of fluid composition information:"
                     f"{', '.join([c.label for c in all_connections])}."
                 )
@@ -981,6 +986,8 @@ def _propagate_fluid_wrappers(self):
                     c.fluid.engine[f] = engine
                 for f, back_end in back_ends.items():
                     c.fluid.back_end[f] = back_end
+                for f, w_kwargs in wrapper_kwargs.items():
+                    c.fluid.wrapper_kwargs[f] = w_kwargs
 
                 c._create_fluid_wrapper()
 

src/tespy/tools/data_containers.py

@@ -802,6 +802,7 @@ def attr():
             "wrapper": dict(),
             "back_end": dict(),
             "engine": dict(),
+            "wrapper_kwargs": dict(),
             "description": None,
             "quantity": None,
             "_is_var": set(),

src/tespy/tools/fluid_properties/__init__.py

@@ -4,6 +4,7 @@
 from .functions import T_mix_ph  # noqa: F401
 from .functions import T_mix_ps  # noqa: F401
 from .functions import T_sat_p  # noqa: F401
+from .functions import conductivity_mix_ph  # noqa: F401
 from .functions import dh_mix_dpQ  # noqa: F401
 from .functions import dT_mix_dph  # noqa: F401
 from .functions import dT_mix_pdh  # noqa: F401
@@ -22,3 +23,4 @@
 from .functions import viscosity_mix_pT  # noqa: F401
 from .helpers import single_fluid  # noqa: F401
 from .wrappers import CoolPropWrapper  # noqa: F401
+from .wrappers import IncompressibleFluidWrapper  # noqa: F401

src/tespy/tools/fluid_properties/functions.py

@@ -352,3 +352,15 @@ def viscosity_mix_pT(p, T, fluid_data, mixing_rule=None):
     else:
         _check_mixing_rule(mixing_rule, V_MIX_PT_DIRECT, "viscosity")
         return VISCOSITY_MIX_PT_DIRECT[mixing_rule](p, T, fluid_data)
+
+
+def conductivity_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None):
+    if get_number_of_fluids(fluid_data) == 1:
+        pure_fluid = get_pure_fluid(fluid_data)
+        return pure_fluid["wrapper"].conductivity_ph(p, h)
+    else:
+        msg = (
+            "Calculation of thermal conductivity is not implemented for "
+            "TESPy based mixtures. You are happily invited to contribute it!"
+        )
+        raise NotImplementedError(msg)

src/tespy/tools/fluid_properties/helpers.py

@@ -360,3 +360,37 @@ def colebrook(reynolds, ks, diameter, darcy_friction_factor, **kwargs):
             / (3.71 * diameter)
         ) + 1 / darcy_friction_factor ** 0.5
     )
+
+
+def _check_fitting_data_structure(x: np.ndarray, y: np.ndarray) -> None:
+    if len(x) != len(y):
+        msg = ""
+        raise ValueError(msg)
+    elif len(x) < 2:
+        msg = ""
+        raise ValueError(msg)
+
+
+def fit_incompressible_viscosity(temperature: np.ndarray, viscosity: np.ndarray) -> tuple:
+    _check_fitting_data_structure(temperature, viscosity)
+
+    x = 1.0 / temperature
+    y = np.log(viscosity)
+
+    return np.polyfit(x, y, 3)
+
+
+def _fit_arrhenius(temperature: np.ndarray, viscosity: np.ndarray) -> tuple:
+    _check_fitting_data_structure(temperature, viscosity)
+    x = 1.0 / temperature
+    y = np.log(viscosity)
+
+    intercept, slope = np.polyfit(x, y, 1)
+
+    return intercept, np.exp(slope)
+
+
+def fit_incompressible_linear(temperature: np.ndarray, y: np.ndarray) -> tuple:
+    _check_fitting_data_structure(temperature, y)
+
+    return np.polyfit(temperature, y, 1)