Fluid properties

The default fluid property engine CoolProp. All available fluids can be found on their homepage. Also see [11]. 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 IncompressibleFluidWrapper. See this section 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 library and

• the 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:

• HEOS: Helmhotz Equation Of State with highest accuracy and lowest calculation speed. This is the default back end!

• REFPROP: Highest accuracy and highest convergence stability. This back end is not free, a separate REFPROP license is required.

• BICUBIC: Tabular back end with high accuracy and very high calculation speed.

• TTSE: Tabular back end with lowest accuracy and very high calculation speed.

• INCOMP: Back end for incompressible fluids.

• IF97: Back end for the IAPWS-IF97 of water, very accurate and much higher calculation speed than HEOS. Due to a bug in the CoolProp back end this option is instable, for more information see the CoolProp issue #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. 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 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 IncompressibleFluidWrapper. See this section for more information.

Fluid mixtures

TESPy provides support for the following built-in mixture rules:

• ideal-cond: gaseous mixtures with flash calculations for water (default).

• ideal: gaseous mixtures without flash calculations.

• incompressible: mass-weighted mixtures of CoolProp incompressible fluids.

• humidair: humid-air mixtures via CoolProp’s HAPropsSI.

• forced-gas: like 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 here. 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.

IncompressibleFluidWrapper

You can use the 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 IncompressibleFluidWrapper will automatically fit functions to your data:

• density and heat capacity through linear interpolation: \(f\left(T\right) = A + B \cdot T\)

• viscosity through an exponential polynomial equation \(\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.

>>> 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.

>>> 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 temperature_data, density_data, heat_capacity_data and viscosity_data!

Then, we can specify the fluid on one of our connections by providing two additional keywords:

• fluid_engines indicating which fluid uses which wrapper class.

• fluid_wrapper_kwargs providing for each fluid (if applicable) the required data.

>>> 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.

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

>>> 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 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 FluidPropertyWrapper class. Below we will use the polynomial formulation for gaseous water from [12] as an example. First we import the necessary dependencies.

>>> 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 kwargs from the connection specification to your concrete wrapper instance, see the section on incompressible fluids.

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

>>> 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.

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

>>> 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 in TESPy

Mixing rules define how pure-fluid properties are combined into mixture properties. You select one per subnetwork via the mixing_rule argument. ideal-cond is the default. The built-in rules are listed in the Fluid Mixtures section above; the underlying equations are documented in the fluid_properties.mixtures module.

Custom mixing rules

New mixing rules can be registered at runtime through the MIXING_RULES registry. Each property function must accept (p, T, fluid_data, **kwargs) and return a scalar value in SI units.

>>> 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 T_ph_inversion and T_ps_inversion keyword arguments (both default to True) control whether the registered h_pT and s_pT functions are also used as Newton residuals for the T(p, h) and T(p, s) inversions.