diff --git a/docs/advanced_features/fluid_properties.rst b/docs/advanced_features/fluid_properties.rst index 228412431..a947f03d1 100644 --- a/docs/advanced_features/fluid_properties.rst +++ b/docs/advanced_features/fluid_properties.rst @@ -4,8 +4,13 @@ Fluid properties ================ The default fluid property engine `CoolProp `_. 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 ` 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 `_ library and - the `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 ` 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 `. + +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 `. .. 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, diff --git a/docs/knowledge_center/faq.rst b/docs/knowledge_center/faq.rst index 47746ba94..37e34fe07 100644 --- a/docs/knowledge_center/faq.rst +++ b/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 ` 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 `. 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 diff --git a/docs/whats_new/v0-9-12.rst b/docs/whats_new/v0-9-12.rst index 44a9d5969..35a7badbe 100644 --- a/docs/whats_new/v0-9-12.rst +++ b/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 ` + (`PR #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 ` + + (`PR #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 `__). - Contributors ############ - Francesco Witte (`@fwitte `__) diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index e48962192..afb858b10 100644 --- a/src/tespy/connections/connection.py +++ b/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): """ diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index ed4ae19d5..e1b09f6b9 100644 --- a/src/tespy/networks/network.py +++ b/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() diff --git a/src/tespy/tools/data_containers.py b/src/tespy/tools/data_containers.py index a63d5615d..778be3ed9 100644 --- a/src/tespy/tools/data_containers.py +++ b/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(), diff --git a/src/tespy/tools/fluid_properties/__init__.py b/src/tespy/tools/fluid_properties/__init__.py index 6351f42e7..dc95484b7 100644 --- a/src/tespy/tools/fluid_properties/__init__.py +++ b/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 diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py index 835bd9d20..00a0bf3c4 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/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) diff --git a/src/tespy/tools/fluid_properties/helpers.py b/src/tespy/tools/fluid_properties/helpers.py index 033ed1801..77dc3700e 100644 --- a/src/tespy/tools/fluid_properties/helpers.py +++ b/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) diff --git a/src/tespy/tools/fluid_properties/wrappers.py b/src/tespy/tools/fluid_properties/wrappers.py index 083f9cd1a..1c2a7f80c 100644 --- a/src/tespy/tools/fluid_properties/wrappers.py +++ b/src/tespy/tools/fluid_properties/wrappers.py @@ -10,9 +10,12 @@ 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 @@ -37,7 +40,7 @@ def __reduce__(self): @wrapper_registry class FluidPropertyWrapper: - def __init__(self, fluid, back_end=None) -> None: + def __init__(self, fluid, back_end=None, **kwargs) -> None: """Base class for fluid property wrappers Parameters @@ -71,7 +74,7 @@ def T_ps(self, p, s): def h_pT(self, p, T): self._not_implemented() - def h_ps(self, p, T): + def h_ps(self, p, s): self._not_implemented() def h_QT(self, Q, T): @@ -122,6 +125,12 @@ def viscosity_ph(self, p, h): def viscosity_pT(self, p, T): self._not_implemented() + def conductivity_ph(self, p, h): + self._not_implemented() + + def conductivity_pT(self, p, T): + self._not_implemented() + def s_ph(self, p, h): self._not_implemented() @@ -132,7 +141,7 @@ def s_pT(self, p, T): @wrapper_registry class CoolPropWrapper(FluidPropertyWrapper): - def __init__(self, fluid, back_end=None) -> None: + def __init__(self, fluid, back_end=None, **kwargs) -> None: """Wrapper for CoolProp.CoolProp.AbstractState instance calls Parameters @@ -356,6 +365,14 @@ def viscosity_pT(self, p, T): self.AS.update(CP.PT_INPUTS, p, T) return self.AS.viscosity() + def conductivity_ph(self, p, h): + self.AS.update(CP.HmassP_INPUTS, h, p) + return self.AS.conductivity() + + def conductivity_pT(self, p, T): + self.AS.update(CP.PT_INPUTS, p, T) + return self.AS.conductivity() + def s_ph(self, p, h): self.AS.update(CP.HmassP_INPUTS, h, p) return self.AS.smass() @@ -365,11 +382,283 @@ def s_pT(self, p, T): return self.AS.smass() +@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 + + 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 + + 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"]) + ) + + 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 + ) + + def h_ps(self, p, s): + return self.h_pT(p, self.T_ps(p, s)) + + def s_ph(self, p, h): + return self.s_pT(p, self.T_ph(p, h)) + + 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) + ) + + 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) + + def T_ps(self, p, s): + return brentq( + self._inverse_s_pT, + self._T_min, + self._T_max, + args=(p, s) + ) + + def conductivity_ph(self, p, h): + return self.conductivity_pT(p, self.T_ph(p, h)) + + def conductivity_pT(self, p, T): + return self._conductivity["A"] * T + self._conductivity["B"] + + def d_ph(self, p, h): + return self.d_pT(p, self.T_ph(p, h)) + + def d_pT(self, p, T): + return self._density["A"] * T + self._density["B"] + + def viscosity_ph(self, p, h): + return self.viscosity_pT(p, self.T_ph(p, h)) + + 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"] + ) + + @wrapper_registry class IAPWSWrapper(FluidPropertyWrapper): - def __init__(self, fluid, back_end=None) -> None: + def __init__(self, fluid, back_end=None, **kwargs) -> None: """Wrapper for iapws library calls Parameters @@ -459,8 +748,8 @@ def phase_ph(self, p, h): return "g" elif phase in ["Two phases", "Saturated vapor", "Saturated liquid"]: return "tp" - else: # to ensure consistent behaviour to CoolPropWrapper - return "phase not recognised" + else: # to ensure consistent behavior to CoolPropWrapper + return "phase not recognized" def d_ph(self, p, h): return self.AS(h=h / 1e3, P=p / 1e6).rho @@ -487,7 +776,7 @@ def s_pT(self, p, T): @wrapper_registry class PyromatWrapper(FluidPropertyWrapper): - def __init__(self, fluid, back_end=None) -> None: + def __init__(self, fluid, back_end=None, **kwargs) -> None: """Wrapper for the Pyromat fluid property library Parameters diff --git a/tests/test_connections.py b/tests/test_connections.py index 9fe6c5458..86401678a 100644 --- a/tests/test_connections.py +++ b/tests/test_connections.py @@ -33,6 +33,7 @@ from tespy.tools.data_containers import FluidProperties as dc_prop from tespy.tools.fluid_properties.functions import T_bubble_p from tespy.tools.fluid_properties.functions import T_dew_p +from tespy.tools.fluid_properties.wrappers import FluidPropertyWrapper from tespy.tools.units import SI_UNITS @@ -527,6 +528,7 @@ def make_connection(cls): QUANTITY_EXEMPTIONS = {} + def properties_of(instance): return [ prop @@ -534,6 +536,7 @@ def properties_of(instance): if isinstance(container, dc_prop) ] + def pytest_generate_tests(metafunc): if "cls_name" in metafunc.fixturenames and "prop" in metafunc.fixturenames: params = [] @@ -543,6 +546,7 @@ def pytest_generate_tests(metafunc): params.append(pytest.param(name, prop, id=f"{cls.__name__}::{prop}")) metafunc.parametrize("cls_name,prop", params) + def test_property_value_not_none(cls_name, prop): instance = make_connection(connection_registry.items[cls_name]) @@ -554,3 +558,17 @@ def test_property_value_not_none(cls_name, prop): ) assert condition, f"Quantity for {prop} of {cls_name} must not be None" + + +def test_wrapper_kwargs_injection(): + so = Source("source") + si = Sink("sink") + c = Connection(so, "out1", si, "in1", label="c") + c.set_attr( + fluid={"H2O": 1}, + fluid_wrapper_kwargs={"H2O": {"testkeyword": "data"}}, + fluid_engines={"H2O": FluidPropertyWrapper} + ) + # if kwargs could not be passed to the Wrapper instantiation this would + # raise an error + c._create_fluid_wrapper() diff --git a/tests/test_networks/test_network.py b/tests/test_networks/test_network.py index 7c1ccfe93..1dd620423 100644 --- a/tests/test_networks/test_network.py +++ b/tests/test_networks/test_network.py @@ -35,6 +35,8 @@ from tespy.connections import Ref from tespy.networks import Network from tespy.tools.data_containers import ComponentMandatoryConstraints as dc_cmc +from tespy.tools.fluid_properties import conductivity_mix_ph +from tespy.tools.fluid_properties.wrappers import IncompressibleFluidWrapper from tespy.tools.helpers import TESPyNetworkError from tespy.tools.helpers import _numeric_deriv @@ -1177,3 +1179,45 @@ class FakeSource(Source): nw.add_conns(c1) nw.solve("design", init_only=True) + + +def test_fluid_kwargs_propagation(): + nw = Network() + nw.units.set_defaults(temperature="°C", pressure="bar") + + pipe = SimpleHeatExchanger("pipe") + + so = Source("source") + si = Sink("sink") + + c1 = Connection(so, "out1", pipe, "in1", label="c1") + c2 = Connection(pipe, "out1", si, "in1", label="c2") + + nw.add_conns(c1, c2) + + fluid_kwargs = { + "temperature_data": np.array([273.15, 373.15]), + "density_data": np.array([1000, 1100]), + "heat_capacity_data": np.array([4000, 4100]), + "viscosity_data": np.array([0.05, 0.00025]), + "conductivity_data": np.array([0.1425, 0.135]) + } + + c1.set_attr( + fluid={"f": 1}, + fluid_engines={"f": IncompressibleFluidWrapper}, + fluid_wrapper_kwargs={"f": fluid_kwargs}, + p=1, T=30 + ) + c2.set_attr(p=0.9, T=50) + pipe.set_attr(Q=1500) + + nw.solve("design") + + # 50 °C is exactly half of range + # heat capacity is implicitly tested, as it is required to find the + # temperature from the enthalpy passed into the function + assert approx(1 / c2.calc_vol()) == 1050 + assert approx( + conductivity_mix_ph(c2.p.val_SI, c2.h.val_SI, c2.fluid_data) + ) == 0.13875 diff --git a/tests/test_tools/test_fluid_properties/test_incompressible.py b/tests/test_tools/test_fluid_properties/test_incompressible.py new file mode 100644 index 000000000..06b97e7f4 --- /dev/null +++ b/tests/test_tools/test_fluid_properties/test_incompressible.py @@ -0,0 +1,135 @@ +# -*- coding: utf-8 + +"""Module for testing fluid properties of incompressibles. + +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 +tests/test_tools/test_fluid_properties/test_incompressible.py + +SPDX-License-Identifier: MIT +""" +import numpy as np +from pytest import approx +from pytest import fixture + +from tespy.tools.fluid_properties.wrappers import IncompressibleFluidWrapper + + +@fixture +def property_data(): + # sample data from a publicly available datasheet + # https://petrocanadalubricants.com/api/sitecore/lubesapi/downloadresource?docID=IM-7852E&type=TechData&lang=english&name=CALFLO%20AF + return { + "temperature_data": np.array([292.647, 310.808, 366.241, 421.673, 477.108, 532.542, 588.826, 618.580]), + "heat_capacity_data": np.array([1901.775, 1961.529, 2143.908, 2326.287, 2508.674, 2691.060, 2876.242, 2974.135]) * 1000, + "density_data": np.array([863.811, 852.596, 818.368, 784.139, 749.909, 715.678, 680.924, 662.551]), + "viscosity_data": np.array([0.050335, 0.028525, 0.007075, 0.002500, 0.00111, 0.000579, 0.000334, 0.000259]), + "conductivity_data": np.array([0.1419, 0.1410, 0.1382, 0.1354 , 0.1327, 0.1299, 0.1271, 0.1256]) + } + + +@fixture +def wrapper_instance(property_data): + return IncompressibleFluidWrapper( + "fluid name", + None, + **property_data + ) + + +def test_setup_wrapper(wrapper_instance): + assert wrapper_instance is not None + + +def test_enthalpy_forwards_backwards(property_data, wrapper_instance): + temperature_data = property_data["temperature_data"] + temperature_check = [] + + for temperature in temperature_data: + h = wrapper_instance.h_pT(None, temperature) + temperature_check.append(wrapper_instance.T_ph(None, h)) + + np.testing.assert_allclose(temperature_check, temperature_data) + + +def test_entropy_forwards_backwards(property_data, wrapper_instance): + temperature_data = property_data["temperature_data"] + temperature_check = [] + + for temperature in temperature_data: + s = wrapper_instance.s_pT(1e5, temperature) + temperature_check.append(wrapper_instance.T_ps(1e5, s)) + + np.testing.assert_allclose(temperature_check, temperature_data) + + +def test_entropy_enthalpy_roundtrip(property_data, wrapper_instance): + temperature_data = property_data["temperature_data"] + temperature_check = [] + + for temperature in temperature_data: + s = wrapper_instance.s_pT(1e5, temperature) + h = wrapper_instance.h_ps(1e5, s) + temperature_check.append(wrapper_instance.T_ph(1e5, h)) + + np.testing.assert_allclose(temperature_check, temperature_data) + + +def test_isentropic(property_data, wrapper_instance): + temperature_data = property_data["temperature_data"] + + enthalpy_inflow = wrapper_instance.h_pT(None, temperature_data) + enthalpy_outflow = wrapper_instance.isentropic(1e5, enthalpy_inflow, 2e5) + rho = wrapper_instance.d_pT(None, temperature_data) + + np.testing.assert_allclose( + enthalpy_outflow - enthalpy_inflow, 1e5 / rho + ) + + +def test_density(property_data, wrapper_instance): + + density_data = property_data["density_data"] + temperature_data = property_data["temperature_data"] + + density = wrapper_instance.d_pT(None, temperature_data) + np.testing.assert_allclose(density, density_data, rtol=1e-3) + + +def test_heat_capacity(property_data, wrapper_instance): + + capacity_data = property_data["heat_capacity_data"] + temperature_data = property_data["temperature_data"] + # capacity is not implemented as method in wrapper, using enthalpy over + # temperature differences + d = 1e-3 + capacity = [ + (h_u - h_l) / (2 * d) + for h_l, h_u in + zip( + wrapper_instance.h_pT(None, temperature_data - d), + wrapper_instance.h_pT(None, temperature_data + d) + ) + ] + np.testing.assert_allclose(capacity, capacity_data, rtol=1e-3) + + +def test_viscosity(property_data, wrapper_instance): + + viscosity_data = property_data["viscosity_data"] + temperature_data = property_data["temperature_data"] + + viscosity = wrapper_instance.viscosity_pT(None, temperature_data) + # allow higher tolerance for viscosity + np.testing.assert_allclose(viscosity, viscosity_data, rtol=1e-2) + + +def test_conductivity(property_data, wrapper_instance): + + conductivity_data = property_data["conductivity_data"] + temperature_data = property_data["temperature_data"] + + conductivity = wrapper_instance.conductivity_pT(None, temperature_data) + # allow higher tolerance for viscosity + np.testing.assert_allclose(conductivity, conductivity_data, rtol=1e-3)