Source code for tespy.tools.fluid_properties.wrappers

# -*- coding: utf-8

"""Module for fluid property wrappers.


This file is part of project TESPy (github.com/oemof/tespy). It's copyrighted
by the contributors recorded in the version control history of the file,
available from its original location
tespy/tools/fluid_properties/wrappers.py

SPDX-License-Identifier: MIT
"""
import CoolProp as CP
import numpy as np
from scipy.optimize import brentq

from tespy.tools.fluid_properties.helpers import fit_incompressible_linear
from tespy.tools.fluid_properties.helpers import fit_incompressible_viscosity
from tespy.tools.global_vars import ERR


[docs] def wrapper_registry(type): wrapper_registry.items[type.__name__] = type return type
wrapper_registry.items = {}
[docs] class SerializableAbstractState(CP.AbstractState): def __init__(self, back_end, fluid_name): self.back_end = back_end self.fluid_name = fluid_name def __reduce__(self): return (self.__class__, (self.back_end, self.fluid_name))
[docs] @wrapper_registry class FluidPropertyWrapper: def __init__(self, fluid, back_end=None, **kwargs) -> None: """Base class for fluid property wrappers Parameters ---------- fluid : str Name of the fluid. back_end : str, optional Name of the back end, by default None """ self.back_end = back_end self.fluid = fluid self.mixture_type = None def _not_implemented(self) -> None: raise NotImplementedError( f"Method is not implemented for {self.__class__.__name__}." )
[docs] def isentropic(self, p_1, h_1, p_2): self._not_implemented()
def _is_below_T_critical(self, T): self._not_implemented()
[docs] def T_ph(self, p, h): self._not_implemented()
[docs] def T_ps(self, p, s): self._not_implemented()
[docs] def h_pT(self, p, T): self._not_implemented()
[docs] def h_ps(self, p, s): self._not_implemented()
[docs] def h_QT(self, Q, T): self._not_implemented()
[docs] def h_pQ(self, p, Q): self._not_implemented()
[docs] def s_QT(self, Q, T): self._not_implemented()
[docs] def T_sat(self, p): self._not_implemented()
[docs] def T_dew(self, p): return self.T_sat(p)
[docs] def T_bubble(self, p): return self.T_sat(p)
[docs] def p_sat(self, T): self._not_implemented()
[docs] def p_dew(self, T): return self.p_sat(T)
[docs] def p_bubble(self, T): return self.p_sat(T)
[docs] def Q_ph(self, p, h): self._not_implemented()
[docs] def phase_ph(self, p, h): self._not_implemented()
[docs] def d_ph(self, p, h): self._not_implemented()
[docs] def d_pT(self, p, T): self._not_implemented()
[docs] def d_QT(self, Q, T): self._not_implemented()
[docs] def viscosity_ph(self, p, h): self._not_implemented()
[docs] def viscosity_pT(self, p, T): self._not_implemented()
[docs] def conductivity_ph(self, p, h): self._not_implemented()
[docs] def conductivity_pT(self, p, T): self._not_implemented()
[docs] def s_ph(self, p, h): self._not_implemented()
[docs] def s_pT(self, p, T): self._not_implemented()
[docs] @wrapper_registry class CoolPropWrapper(FluidPropertyWrapper): def __init__(self, fluid, back_end=None, **kwargs) -> None: """Wrapper for CoolProp.CoolProp.AbstractState instance calls Parameters ---------- fluid : str Name of the fluid back_end : str, optional CoolProp back end for the AbstractState object, by default "HEOS" """ super().__init__(fluid, back_end) if self.back_end is None: self.back_end = "HEOS" self._identify_mixture() self.AS = SerializableAbstractState(self.back_end, self.fluid) self._set_mixture_fractions() self._set_constants() self._last_ip = None self._last_a = None self._last_b = None def _update(self, input_pair, a, b): if input_pair == self._last_ip and a == self._last_a and b == self._last_b: return self._last_ip = None self.AS.update(input_pair, a, b) self._last_ip = input_pair self._last_a = a self._last_b = b def _identify_mixture(self): """Parse the fluid name to identify, if and what kind of mixture we are working with """ if "[" in self.fluid: if "|" not in self.fluid: msg = ( f"The fluid {self.fluid} requires the specification of " "mass, volume or molar based composition information." "You can do this by appending '|' and 'mass' at the end " "of the fluid string. For example, " "'NAMEOFFLUID[0.5]|mass' to indicate a mass based mixture." ) raise ValueError(msg) self.fluid, self.mixture_type = self.fluid.split("|") allowed = ["mass", "molar", "volume"] if self.mixture_type not in allowed: msg = ( "For the specification of the composition type you have " f"to select from {', '.join(allowed)}." ) if "&" in self.fluid: _fluids_with_fractions = self.fluid.split("&") else: _fluids_with_fractions = [self.fluid] fluid_names = [] fractions = [] for fluid in _fluids_with_fractions: if "[" in fluid: _fluid_name, _fraction = fluid.split("[") _fraction = float(_fraction.replace("]", "")) fractions += [_fraction] else: _fluid_name = fluid fluid_names += [_fluid_name] self.fractions = fractions self.fluid = "&".join(fluid_names) def _set_mixture_fractions(self): """Set the fractions for provided mixture""" if self.mixture_type == "mass": self.AS.set_mass_fractions(self.fractions) elif self.mixture_type == "molar": self.AS.set_mole_fractions(self.fractions) elif self.mixture_type == "volume": self.AS.set_volu_fractions(self.fractions) def _set_constants(self): """Setup constants for later quick access, e.g. mixture fractions minimum/maximum pressure/temperature and critical point properties """ self._T_min = self.AS.trivial_keyed_output(CP.iT_min) self._T_max = self.AS.trivial_keyed_output(CP.iT_max) if self.back_end == "INCOMP": self._p_min = 1e2 self._p_max = 1e8 self._p_crit = 1e8 self._T_crit = None self._molar_mass = 1 if self.mixture_type is not None: try: self._T_min = max( self.AS.trivial_keyed_output(CP.iT_freeze), self._T_min ) except ValueError: pass else: if self.back_end == "HEOS": # see https://github.com/CoolProp/CoolProp/discussions/2443 self._T_max *= 1.45 if self.back_end == "REFPROP": if self.mixture_type is not None: self._T_min += 5 self._p_min = 1e1 else: self._p_min = self.AS.trivial_keyed_output(CP.iP_min) self._p_max = self.AS.trivial_keyed_output(CP.iP_max) self._p_crit = self.AS.trivial_keyed_output(CP.iP_critical) self._T_crit = self.AS.trivial_keyed_output(CP.iT_critical) self._molar_mass = self.AS.trivial_keyed_output(CP.imolar_mass) def _is_below_T_critical(self, T): return T < self._T_crit
[docs] def get_T_max(self, p): if self.back_end == "INCOMP": return self.T_sat(p) else: return self._T_max
[docs] def isentropic(self, p_1, h_1, p_2): return self.h_ps(p_2, self.s_ph(p_1, h_1))
[docs] def T_ph(self, p, h): self._update(CP.HmassP_INPUTS, h, p) return self.AS.T()
[docs] def T_ps(self, p, s): self._update(CP.PSmass_INPUTS, p, s) return self.AS.T()
[docs] def h_pQ(self, p, Q): self._update(CP.PQ_INPUTS, p, Q) return self.AS.hmass()
[docs] def h_ps(self, p, s): self._update(CP.PSmass_INPUTS, p, s) return self.AS.hmass()
[docs] def h_pT(self, p, T): self._update(CP.PT_INPUTS, p, T) return self.AS.hmass()
[docs] def h_QT(self, Q, T): self._update(CP.QT_INPUTS, Q, T) return self.AS.hmass()
[docs] def s_QT(self, Q, T): self._update(CP.QT_INPUTS, Q, T) return self.AS.smass()
[docs] def T_sat(self, p): self._update(CP.PQ_INPUTS, p, 0) return self.AS.T()
[docs] def T_dew(self, p): self._update(CP.PQ_INPUTS, p, 1) return self.AS.T()
[docs] def T_bubble(self, p): self._update(CP.PQ_INPUTS, p, 0) return self.AS.T()
[docs] def p_sat(self, T): self._update(CP.QT_INPUTS, 0.5, T) return self.AS.p()
[docs] def p_dew(self, T): self._update(CP.QT_INPUTS, 1, T) return self.AS.p()
[docs] def p_bubble(self, T): self._update(CP.QT_INPUTS, 0, T) return self.AS.p()
[docs] def Q_ph(self, p, h): self._update(CP.HmassP_INPUTS, h, p) if len(self.fractions) > 1: return self.AS.Q() phase = self.AS.phase() if phase == CP.iphase_twophase: return self.AS.Q() elif phase == CP.iphase_liquid: return 0 elif phase == CP.iphase_gas: return 1 else: # all other phases - though this should be unreachable as p is sub-critical return -1
[docs] def phase_ph(self, p, h): if self.back_end == "INCOMP": return "state not recognized" self._update(CP.HmassP_INPUTS, h, p) phase = self.AS.phase() if phase == CP.iphase_twophase: return "tp" elif phase == CP.iphase_liquid: return "l" elif phase == CP.iphase_gas: return "g" elif phase == CP.iphase_supercritical_gas: return "g" else: return "state not recognised"
[docs] def d_ph(self, p, h): self._update(CP.HmassP_INPUTS, h, p) return self.AS.rhomass()
[docs] def d_pT(self, p, T): self._update(CP.PT_INPUTS, p, T) return self.AS.rhomass()
[docs] def d_QT(self, Q, T): self._update(CP.QT_INPUTS, Q, T) return self.AS.rhomass()
[docs] def viscosity_ph(self, p, h): self._update(CP.HmassP_INPUTS, h, p) return self.AS.viscosity()
[docs] def viscosity_pT(self, p, T): self._update(CP.PT_INPUTS, p, T) return self.AS.viscosity()
[docs] def conductivity_ph(self, p, h): self._update(CP.HmassP_INPUTS, h, p) return self.AS.conductivity()
[docs] def conductivity_pT(self, p, T): self._update(CP.PT_INPUTS, p, T) return self.AS.conductivity()
[docs] def s_ph(self, p, h): self._update(CP.HmassP_INPUTS, h, p) return self.AS.smass()
[docs] def s_pT(self, p, T): self._update(CP.PT_INPUTS, p, T) return self.AS.smass()
[docs] @wrapper_registry class IncompressibleFluidWrapper(FluidPropertyWrapper): """Class to represent a fluid in TESPy using tabular data Parameters ---------- fluid : str Name of fluid back_end : str, optional Name of the back end in context of CoolProp, by default None temperature_data : np.ndarray Array of temperature measurements in SI units (Kelvin) density_data : np.ndarray Array of corresponding density values in SI units (kg/m3) heat_capacity_data : np.ndarray Array of corresponding heat capacity values in SI units (J/kg) viscosity_data : np.ndarray Array of corresponding **dynamic** viscosity values in SI units (Pas) conductivity_data : np.ndarray Array of corresponding thermal conductivity values in SI units (W/mK) """ def __init__(self, fluid, back_end=None, **kwargs): """Class to represent a fluid in TESPy using tabular data Parameters ---------- fluid : str Name of fluid back_end : str, optional Name of the back end in context of CoolProp, by default None temperature_data : np.ndarray Array of temperature measurements in SI units (Kelvin) density_data : np.ndarray Array of corresponding density values in SI units (kg/m3) heat_capacity_data : np.ndarray Array of corresponding heat capacity values in SI units (J/kg) viscosity_data : np.ndarray Array of corresponding **dynamic** viscosity values in SI units (Pas) conductivity_data : np.ndarray Array of corresponding thermal conductivity values in SI units (W/mK) """ super().__init__(fluid, back_end, **kwargs) self.temperature_data = None self.heat_capacity_data = None self.density_data = None self.viscosity_data = None for key in ["temperature", "heat_capacity", "density", "viscosity"]: value = kwargs.get(f"{key}_data") if value is None: msg = ( f"The {self.__class__.__name__} requires specification of " f"the '{key}_data' keyword in the form of a numpy array." ) raise KeyError(msg) else: setattr(self, f"{key}_data", value) self.conductivity_data = kwargs.get("conductivity_data") self._T_ref = kwargs.get("T_ref", min(self.temperature_data)) self._p_ref = kwargs.get("p_ref", 1e5) self._fit_data() self._set_constants() def _fit_data(self): A, B = fit_incompressible_linear( self.temperature_data, self.heat_capacity_data ) self._heat_capacity = { "A": A, "B": B } A, B = fit_incompressible_linear( self.temperature_data, self.density_data ) self._density = { "A": A, "B": B } if self.conductivity_data is not None: A, B = fit_incompressible_linear( self.temperature_data, self.conductivity_data ) else: A, B = np.nan, np.nan self._conductivity = { "A": A, "B": B } A, B, C, D = fit_incompressible_viscosity( self.temperature_data, self.viscosity_data ) self._viscosity = { "A": A, "B": B, "C": C, "D": D } def _set_constants(self): # evaluate h at T=T_ref self._h_ref = self._h_pT(None, self._T_ref) self._T_min = self._T_ref self._T_max = max(self.temperature_data) self._molar_mass = 1 self._p_min = 100 self._p_max = 10000000 self._p_crit = self._p_max self._T_crit = None
[docs] def get_fitting_report(self): import matplotlib.pyplot as plt def plot_property(ax, temperature, measurements, evaluation): _fit, = ax.plot(temperature, evaluation, "-", color="red") _data = ax.scatter(temperature, measurements, marker="x", c="blue") ax_err = ax.twinx() _err = ax_err.scatter( temperature, (evaluation - measurements) / measurements * 100, c="#0000ff66" ) ax_err.set_ylabel("Deviation between fit and data in %") return [_data, _fit, _err] fig, ax = plt.subplots(2, 2, figsize=(10, 10), sharex=True) ax[0, 0].set_title("Heat capacity") ax[0, 1].set_title("Density") ax[1, 0].set_title("Viscosity") ax[1, 1].set_title("Thermal conductivity") temperature_data = self.temperature_data heat_capacity_data = self.heat_capacity_data density_data = self.density_data viscosity_data = self.viscosity_data conductivity_data = self.conductivity_data d = 0.001 heat_capacity_eval = ( self.h_pT(None, temperature_data + d) - self.h_pT(None, temperature_data - d) ) / (2 * d) density_eval = self.d_pT(None, temperature_data) viscosity_eval = self.viscosity_pT(None, temperature_data) conductivity_eval = self.conductivity_pT(None, temperature_data) lines = plot_property(ax[0, 0], temperature_data, heat_capacity_data, heat_capacity_eval) labels = ["datapoints", "fitted function", "deviation"] plot_property(ax[0, 1], temperature_data, density_data, density_eval) plot_property(ax[1, 0], temperature_data, viscosity_data, viscosity_eval) ax[1, 0].set_yscale("log") plot_property(ax[1, 1], temperature_data, conductivity_data, conductivity_eval) ax[0, 0].set_ylabel("Heat capacity in J/kgK") ax[0, 1].set_ylabel("Density in kg/m3") ax[1, 0].set_ylabel("Viscosity in Pas") ax[1, 1].set_ylabel("Thermal conductivity in W/mK") ax[1, 0].set_xlabel("Temperature in K") ax[1, 1].set_xlabel("Temperature in K") fig.legend( lines, labels, loc="upper center", ncol=3, bbox_to_anchor=(0.5, 1.05) ) plt.tight_layout() return fig, ax
[docs] def T_ph(self, p, h): # Inverse function of h_pT, using quadratic formula with adding the # root return ( ( -self._heat_capacity["B"] + ( self._heat_capacity["B"] ** 2 - 4 * 0.5 * self._heat_capacity["A"] * - (h + self._h_ref) ) ** 0.5 ) / (2 * 0.5 * self._heat_capacity["A"]) )
[docs] def h_pT(self, p, T): return self._h_pT(p, T) - self._h_ref
def _h_pT(self, p, T): # h = integral cp(T) dT return ( 0.5 * self._heat_capacity["A"] * T ** 2 + self._heat_capacity["B"] * T )
[docs] def h_ps(self, p, s): return self.h_pT(p, self.T_ps(p, s))
[docs] def s_ph(self, p, h): return self.s_pT(p, self.T_ph(p, h))
[docs] def s_pT(self, p, T): # s0 = 0 return ( self._heat_capacity["B"] * np.log(T / self._T_ref) + self._heat_capacity["A"] * (T - self._T_ref) - self.d_pT(p, T) * (p - self._p_ref) )
[docs] def isentropic(self, p_1, h_1, p_2): # assumption that temperature barely changes T = self.T_ph(p_1, h_1) return h_1 + (p_2 - p_1) / self.d_pT(p_1, T)
def _inverse_s_pT(self, T, p, s): return s - self.s_pT(p, T)
[docs] def T_ps(self, p, s): return brentq( self._inverse_s_pT, self._T_min, self._T_max, args=(p, s) )
[docs] def conductivity_ph(self, p, h): return self.conductivity_pT(p, self.T_ph(p, h))
[docs] def conductivity_pT(self, p, T): return self._conductivity["A"] * T + self._conductivity["B"]
[docs] def d_ph(self, p, h): return self.d_pT(p, self.T_ph(p, h))
[docs] def d_pT(self, p, T): return self._density["A"] * T + self._density["B"]
[docs] def viscosity_ph(self, p, h): return self.viscosity_pT(p, self.T_ph(p, h))
[docs] def viscosity_pT(self, p, T): return np.exp( self._viscosity["A"] / T ** 3 + self._viscosity["B"] / T ** 2 + self._viscosity["C"] / T + self._viscosity["D"] )
[docs] @wrapper_registry class IAPWSWrapper(FluidPropertyWrapper): def __init__(self, fluid, back_end=None, **kwargs) -> None: """Wrapper for iapws library calls Parameters ---------- fluid : str Name of the fluid back_end : str, optional CoolProp back end for the AbstractState object, by default "IF97" """ # avoid unnecessary loading time if not used try: import iapws except ModuleNotFoundError: msg = ( "To use the iapws fluid properties you need to install " "iapws." ) raise ModuleNotFoundError(msg) if back_end is None: back_end = "IF97" super().__init__(fluid, back_end) if self.back_end == "IF97": self.AS = iapws.IAPWS97 elif self.back_end == "IF95": self.AS = iapws.IAPWS95 else: msg = f"The specified back_end {self.back_end} is not available." raise NotImplementedError(msg) self._set_constants(iapws) def _set_constants(self, iapws): self._T_min = iapws._iapws.Tt self._T_max = 2000 self._p_min = iapws._iapws.Pt * 1e6 self._p_max = 100e6 self._p_crit = iapws._iapws.Pc * 1e6 self._T_crit = iapws._iapws.Tc self._molar_mass = iapws._iapws.M def _is_below_T_critical(self, T): return T < self._T_crit
[docs] def isentropic(self, p_1, h_1, p_2): return self.h_ps(p_2, self.s_ph(p_1, h_1))
[docs] def T_ph(self, p, h): return self.AS(h=h / 1e3, P=p / 1e6).T
[docs] def T_ps(self, p, s): return self.AS(s=s / 1e3, P=p / 1e6).T
[docs] def h_pQ(self, p, Q): return self.AS(P=p / 1e6, x=Q).h * 1e3
[docs] def h_ps(self, p, s): return self.AS(P=p / 1e6, s=s / 1e3).h * 1e3
[docs] def h_pT(self, p, T): return self.AS(P=p / 1e6, T=T).h * 1e3
[docs] def h_QT(self, Q, T): return self.AS(T=T, x=Q).h * 1e3
[docs] def s_QT(self, Q, T): return self.AS(T=T, x=Q).s * 1e3
[docs] def T_sat(self, p): return self.AS(P=p / 1e6, x=0).T
[docs] def p_sat(self, T): if T > self._T_crit: T = self._T_crit * 0.99 return self.AS(T=T / 1e6, x=0).P * 1e6
[docs] def Q_ph(self, p, h): return self.AS(h=h / 1e3, P=p / 1e6).x
[docs] def phase_ph(self, p, h): phase = self.AS(h=h / 1e3, P=p / 1e6).phase if phase in ["Liquid"]: return "l" elif phase in ["Vapour"]: return "g" elif phase in ["Two phases", "Saturated vapor", "Saturated liquid"]: return "tp" else: # to ensure consistent behavior to CoolPropWrapper return "phase not recognized"
[docs] def d_ph(self, p, h): return self.AS(h=h / 1e3, P=p / 1e6).rho
[docs] def d_pT(self, p, T): return self.AS(T=T, P=p / 1e6).rho
[docs] def d_QT(self, Q, T): return self.AS(T=T, x=Q).rho
[docs] def viscosity_ph(self, p, h): return self.AS(P=p / 1e6, h=h / 1e3).mu
[docs] def viscosity_pT(self, p, T): return self.AS(T=T, P=p / 1e6).mu
[docs] def s_ph(self, p, h): return self.AS(P=p / 1e6, h=h / 1e3).s * 1e3
[docs] def s_pT(self, p, T): return self.AS(P=p / 1e6, T=T).s * 1e3
[docs] @wrapper_registry class PyromatWrapper(FluidPropertyWrapper): def __init__(self, fluid, back_end=None, **kwargs) -> None: """Wrapper for the Pyromat fluid property library Parameters ---------- fluid : str Name of the fluid back_end : str, optional CoolProp back end for the AbstractState object, by default None """ # avoid unnecessary loading time if not used try: import pyromat as pm pm.config['unit_energy'] = "J" pm.config['unit_pressure'] = "Pa" pm.config['unit_molar'] = "mol" except ModuleNotFoundError: msg = ( "To use the pyromat fluid properties you need to install " "pyromat." ) raise ModuleNotFoundError(msg) super().__init__(fluid, back_end) self._create_AS(pm) self._set_constants() def _create_AS(self, pm): self.AS = pm.get(f"{self.back_end}.{self.fluid}") def _set_constants(self): self._p_min, self._p_max = 100, 1000e5 self._T_min, self._T_max = self.AS.Tlim() if self.back_end == "mp": self._T_crit, self._p_crit = self.AS.critical() else: self._T_crit = self._T_max self._p_crit = self._p_max self._molar_mass = self.AS.mw()
[docs] def isentropic(self, p_1, h_1, p_2): return self.h_ps(p_2, self.s_ph(p_1, h_1))
[docs] def T_ph(self, p, h): return self.AS.T(p=p, h=h)[0]
[docs] def T_ps(self, p, s): return self.AS.T(p=p, s=s)[0]
[docs] def h_pT(self, p, T): return self.AS.h(p=p, T=T)[0]
[docs] def h_ps(self, p, s): return self.AS.h(p=p, s=s)[0]
[docs] def d_ph(self, p, h): return self.AS.d(p=p, h=h)[0]
[docs] def d_pT(self, p, T): return self.AS.d(p=p, T=T)[0]
[docs] def s_ph(self, p, h): return self.AS.s(p=p, h=h)[0]
[docs] def s_pT(self, p, T): return self.AS.s(p=p, T=T)[0]
[docs] def h_QT(self, Q, T): if self.back_end == "ig": self._not_implemented() return self.AS.h(x=Q, T=T)[0]
[docs] def h_pQ(self, p, Q): if self.back_end == "ig": self._not_implemented() return self.AS.h(p=p, x=Q)[0]
[docs] def s_QT(self, Q, T): if self.back_end == "ig": self._not_implemented() return self.AS.s(x=Q, T=T)[0]
[docs] def Q_ph(self, p, h): if self.back_end == "ig": self._not_implemented() return self.AS.x(p=p, h=h)[0]
[docs] def d_QT(self, Q, T): if self.back_end == "ig": self._not_implemented() return self.AS.d(x=Q, T=T)[0]