.. _fluid_properties_label:

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 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 own one, i.e.

- the `iapws <https://github.com/jjgomera/iapws/>`_ library and
- the `pyromat <https://github.com/chmarti1/PYroMat/>`_ library.

For each fluid you can specify, which library should be used, and you can easily
implement your own engine, for example, if your fluid is not available through
the predefined engines, or you need a different implementation of the fluid
properties of a specific fluid. Further below on this page we have added an
example implementation of the KKH polynomial formulation to illustrate, how you
can integrate your own engine into TESPy.

CoolProp
--------
CoolProp is the default fluid property back end of TESPy. It provides different
back ends for different applications and allows calling mixture functions.

Available back ends
+++++++++++++++++++
CoolProp provides multiple back ends for fluid property calculation. The
back ends vary in calculation speed and calculation accuracy. It is possible
to choose from the following back ends:

- :code:`HEOS`: Helmhotz Equation Of State with highest accuracy and lowest
  calculation speed. **This is the default back end!**
- :code:`REFPROP`: Highest accuracy and highest convergence stability.
  **This back end is not free**, a separate
  `REFPROP <https://www.nist.gov/srd/refprop>`__ license is required.
- :code:`BICUBIC`: Tabular back end with high accuracy and very high
  calculation speed.
- :code:`TTSE`: Tabular back end with lowest accuracy and very high calculation
  speed.
- :code:`INCOMP`: Back end for incompressible fluids.
- :code:`IF97`: Back end for the IAPWS-IF97 of water, very accurate and much
  higher calculation speed than :code:`HEOS`. Due to a bug in the CoolProp
  back end this option is instable, for more information see the
  `CoolProp issue #1918 <https://github.com/CoolProp/CoolProp/issues/1918/>`_.

For more information on the back ends please visit the CoolProp online
documentation.

Pure and pseudo-pure fluids
+++++++++++++++++++++++++++
CoolProp covers the most important fluids such as water, air as a pseudo-pure
fluid as well as its components, several fuels and refrigerants etc.. Look for
the aliases in the list of
`fluids <https://coolprop.org/fluid_properties/PurePseudoPure.html>`__.
All fluids provided in this list cover liquid and gaseous state and the
two-phase region.

Incompressible fluids
+++++++++++++++++++++
If you are looking for heat transfer fluids, the list of incompressible
`fluids <https://coolprop.org/fluid_properties/Incompressibles.html>`__
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 the following built-in mixture rules:

- :code:`ideal-cond`: gaseous mixtures **with flash calculations for water** (default).
- :code:`ideal`: gaseous mixtures **without flash calculations**.
- :code:`incompressible`: mass-weighted mixtures of CoolProp incompressible fluids.
- :code:`humidair`: humid-air mixtures via CoolProp's :code:`HAPropsSI`.
- :code:`forced-gas`: like :code:`ideal`, but forces the gas phase for water at saturation.

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>`.  Custom mixing rules can also be
registered at runtime - see the same section for details.

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",
    ...     pressure_difference="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)

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
    ...     "conductivity_data": np.array([0.15, 0.13])  # W/m/K
    ... }

.. 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
-------------------
To use any of the other fluid property engines, you can do the following, e.g.
to use the `iapws` back end:

.. code-block:: python

    >>> from tespy.components import Sink
    >>> from tespy.components import Source
    >>> from tespy.components import Turbine
    >>> from tespy.connections import Connection
    >>> from tespy.networks import Network
    >>> from tespy.tools.fluid_properties.wrappers import IAPWSWrapper

    >>> nwk = Network(iterinfo=False)

    >>> so = Source("Source")
    >>> tu = Turbine("Turbine")
    >>> si = Sink("Sink")

    >>> c1 = Connection(so, "out1", tu, "in1", label="1")
    >>> c2 = Connection(tu, "out1", si, "in1", label="2")

    >>> nwk.add_conns(c1, c2)

    >>> tu.set_attr(eta_s=0.9)

    >>> c1.set_attr(
    ...     v=1, p=1e5, T=500,
    ...     fluid={"IF97::H2O": 1}, fluid_engines={"H2O": IAPWSWrapper}
    ... )
    >>> c2.set_attr(p=1e4)

    >>> nwk.solve("design")
    >>> float(round(c2.x.val, 3))
    0.99

    >>> tu.set_attr(eta_s=None)
    >>> c2.set_attr(x=1)

    >>> nwk.solve("design")
    >>> float(round(tu.eta_s.val, 3))
    0.841


Implementing a custom engine
----------------------------
The fluid property calls to different engines have to be masqueraded with
respective wrappers. The implementation of the wrappers for `CoolProp`,
`iapws` and `pyromat` can be found in the
:py:mod:`fluid_properties.wrappers <tespy.tools.fluid_properties.wrappers>`
module, and serve as example implementations for your own wrappers:

The wrapper for your own engine (or an engine from a different library) has to
inherit from the
:py:class:`FluidPropertyWrapper <tespy.tools.fluid_properties.wrappers.FluidPropertyWrapper>`
class. Below we will use the polynomial formulation for **gaseous water** from
:cite:`Knacke1991` as an example. First we import the necessary dependencies.

.. code-block:: python

    >>> import numpy as np
    >>> from tespy.tools.fluid_properties.wrappers import FluidPropertyWrapper
    >>> from tespy.tools.global_vars import GAS_CONSTANT_UNI

Then we set up a new class and implement the methods to calculate enthalpy and
entropy from (pressure and) temperature. The structure and names of the
functions have to match the pattern from the `FluidPropertyWrapper`, in this
case `h_pT`. On top of that, we add a backwards function `T_ph` and a function
to analytically calculate the heat capacity `_cp_pT`, the derivative of the
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

    # coefficients       H+       S+       a      b       c    d        M
    >>> COEF = {
    ...    "H2O": [-253.871, -11.750, 34.376, 7.841, -0.423, 0.0, 18.0152],
    ... }

    >>> class KKHWrapper(FluidPropertyWrapper):
    ...
    ...     def __init__(self, fluid, back_end=None, reference_temperature=298.15, **kwargs) -> None:
    ...         super().__init__(fluid, back_end)
    ...
    ...         if self.fluid not in COEF:
    ...             msg = "Fluid not available in KKH database"
    ...             raise KeyError(msg)
    ...
    ...         self.coefficients = COEF[fluid]
    ...         self.h_ref = self._h_pT(None, reference_temperature)
    ...         self._molar_mass = self.coefficients[-1] * 1e-3
    ...         self._T_min = 100
    ...         self._T_max = 2000
    ...         self._p_min = 1000
    ...         self._p_max = 10000000
    ...         self._T_crit = 647.096
    ...         self._p_crit = 2206400
    ...
    ...     def cp_pT(self, p, T):
    ...         y = T * 1e-3
    ...         return 1e3 * (
    ...             self.coefficients[2]
    ...             + self.coefficients[3] * y
    ...             + self.coefficients[4] / (y ** 2)
    ...             + self.coefficients[5] * y ** 2
    ...         ) / self.coefficients[6]
    ...
    ...     def h_pT(self, p, T):
    ...         return self._h_pT(p, T) - self.h_ref
    ...
    ...     def _h_pT(self, p, T):
    ...         y = T * 1e-3
    ...         return 1e6 * (
    ...             self.coefficients[0]
    ...             + self.coefficients[2] * y
    ...             + self.coefficients[3] / 2 * y ** 2
    ...             - self.coefficients[4] / y
    ...             + self.coefficients[5] / 3 * y ** 3
    ...         ) / self.coefficients[6]
    ...
    ...     def T_ph(self, p, h):
    ...         return newton(self.h_pT, self.cp_pT, h, p)
    ...
    ...     def isentropic(self, p_1, h_1, p_2):
    ...         T_1 = self.T_ph(p_1, h_1)
    ...         cp = self.cp_pT(p_1, T_1)
    ...         kappa = cp / (cp - GAS_CONSTANT_UNI / self._molar_mass)
    ...         T_2 = T_1 * (p_2 / p_1) ** ((kappa - 1) / kappa)
    ...         return self.h_pT(p_2, T_2)

We can add a newton for the backwards function:

.. code-block:: python

    >>> def newton(func, deriv, h, p):
    ...     # default valaues
    ...     T = 300
    ...     valmin = 70
    ...     valmax = 3000
    ...     max_iter = 10
    ...     tol_rel = 1e-6
    ...
    ...     # start newton loop
    ...     expr = True
    ...     i = 0
    ...     while expr:
    ...         # calculate function residual and new value
    ...         res = h - func(p, T)
    ...         T += res / deriv(p, T)
    ...
    ...         # check for value ranges
    ...         if T < valmin:
    ...             T = valmin
    ...         if T > valmax:
    ...             T = valmax
    ...         i += 1
    ...
    ...         if i > max_iter:
    ...             break
    ...
    ...         expr = abs(res / h) >= tol_rel
    ...
    ...     return T

And then we can test a call to the interface and check our results, also compare
them to CoolProp.

.. code-block:: python

    >>> kkh_water = KKHWrapper("H2O", reference_temperature=298.15)  # same as in CoolProp
    >>> h = kkh_water.h_pT(1e5, 400)
    >>> T = kkh_water.T_ph(1e5, h)
    >>> round(T, 1)
    400.0
    >>> round(h)
    189769

    >>> from tespy.tools.fluid_properties import CoolPropWrapper

    >>> coolprop_water = CoolPropWrapper("H2O")
    >>> h_cp = coolprop_water.h_pT(1e5, 400)
    >>> T_cp = coolprop_water.T_ph(1e5, h_cp)
    >>> round(T_cp, 1)
    400.0
    >>> round(h)
    189769

To use this wrapper in a simple TESPy model, we can then proceed as we have in
the previous section:

.. code-block:: python

    >>> from tespy.components import Sink
    >>> from tespy.components import Source
    >>> from tespy.components import Turbine
    >>> from tespy.connections import Connection
    >>> from tespy.networks import Network

    >>> nwk = Network(iterinfo=False)
    >>> nwk.units.set_defaults(
    ...     temperature="degC", pressure="MPa", pressure_difference="MPa"
    ... )

    >>> so = Source("Source")
    >>> tu = Turbine("Turbine")
    >>> si = Sink("Sink")

    >>> c1 = Connection(so, "out1", tu, "in1", label="1")
    >>> c2 = Connection(tu, "out1", si, "in1", label="2")

    >>> nwk.add_conns(c1, c2)

    >>> c1.set_attr(
    ...     m=1, p=10, T=600,
    ...     fluid={"H2O": 1}, fluid_engines={"H2O": KKHWrapper}
    ... )
    >>> c2.set_attr(p=1, T=400)

    >>> nwk.solve("design")

    >>> tu.set_attr(eta_s=0.9)
    >>> c2.set_attr(T=None)
    >>> nwk.solve("design")
    >>> round(c2.T.val, 1)
    306.3

.. _mixture_routines_label:

Mixture routines in TESPy
-------------------------
Mixing rules define how pure-fluid properties are combined into mixture
properties. You select one per subnetwork via the :code:`mixing_rule` argument.
:code:`ideal-cond` is the default. The built-in rules are listed in the
:ref:`Fluid Mixtures <fluid_properties_label>` section above; the underlying
equations are documented in the
:py:mod:`fluid_properties.mixtures <tespy.tools.fluid_properties.mixtures>`
module.

Custom mixing rules
+++++++++++++++++++
New mixing rules can be registered at runtime through the
:py:obj:`MIXING_RULES <tespy.tools.fluid_properties.MIXING_RULES>` registry.
Each property function must accept :code:`(p, T, fluid_data, **kwargs)` and
return a scalar value in SI units.

.. code-block:: python

    >>> from tespy.tools.fluid_properties import MIXING_RULES

    >>> def my_h_pT(p, T, fluid_data, **kwargs):
    ...     h = 0
    ...     for data in fluid_data.values():
    ...         h += data["wrapper"].h_pT(p, T) * data["mass_fraction"]
    ...     return h

    >>> def my_s_pT(p, T, fluid_data, **kwargs):
    ...     s = 0
    ...     for data in fluid_data.values():
    ...         s += data["wrapper"].s_pT(p, T) * data["mass_fraction"]
    ...     return s

    >>> MIXING_RULES.register(
    ...     "my-rule",
    ...     h_pT=my_h_pT,
    ...     s_pT=my_s_pT,
    ... )

The :code:`T_ph_inversion` and :code:`T_ps_inversion` keyword arguments
(both default to :code:`True`) control whether the registered :code:`h_pT`
and :code:`s_pT` functions are also used as Newton residuals for the
:code:`T(p, h)` and :code:`T(p, s)` inversions.