diff --git a/README.rst b/README.rst index 398ad3096..fefd837ca 100644 --- a/README.rst +++ b/README.rst @@ -124,17 +124,16 @@ use. Get in touch ============ -Online "Stammtisch" -------------------- +Community meetings +------------------ -We have decided to start a reoccurring "Stammtisch" meeting for all interested -TESPy users and (potential) developers. You are invited to join us on every 3rd -Monday of a month at 17:00 CE(S)T for a casual get together. The first meeting -will be held at June, 20, 2022. The intent of this meeting is to establish a -more active and well-connected network of TESPy users and developers. +There are online QA meetings every month, please check the +`oemof calendar `__ for date, time and meeting +link. You are invited to join us with your questions, issues and suggestions! -If you are interested, you can simply join the meeting at -https://meet.jit.si/tespy_user_meeting. We are looking forward to seeing you! +Furthermore, there are in-person community meetings. These are held once or +twice a year. To learn about upcoming meetings, follow the blog over at +``__ for the announcements. User forum ---------- 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/advanced_tutorials/heat_pump_steps.rst b/docs/advanced_tutorials/heat_pump_steps.rst index b3b09ebac..15877507b 100644 --- a/docs/advanced_tutorials/heat_pump_steps.rst +++ b/docs/advanced_tutorials/heat_pump_steps.rst @@ -142,7 +142,7 @@ decides the overall mass flow in the systems. In this tutorial we will first build the system with parameters that ensure stable starting values for a simulation, which in the end will be switched to reasonable values for the individual parts of the system. For - example, instead of the evaporation pressure we will use the terminal + example, instead of the evaporation pressure we will use the pinch temperature difference at the condenser instead. .. literalinclude:: /../tutorial/advanced/stepwise.py @@ -152,11 +152,11 @@ decides the overall mass flow in the systems. In order to calculate this network further parametrization is necessary, as e.g. the fluids are not determined yet: At the hot inlet of the condenser we -define the temperature, pressure and the fluid informaton. A good guess for -pressure can be obtained from CoolProp's PropsSI function. We know that the -condensation temperature must be higher than the consumer's feed flow -temperature. Therefore, we can set the pressure to a slightly higher value of -that temperature's corresponding condensation pressure. +define the temperature and the fluid information. On top we can define the +pressure level by specifying the dew line temperature. We know, that the +condensation temperature will be slightly below the heat demand feed flow +temperature. Therefore, we can set the dew line temperature to a slightly lower +value of that temperature. The same needs to be done for the consumer cycle. We suggest setting the parameters at the pump's inlet. On top, we assume that the consumer requires a @@ -240,7 +240,9 @@ Connections Since the old connection :code:`1` lead to a sink, we have to replace this connection in the network. We can do that by using the method :code:`del_conns` passing :code:`c1`. After that, we can create the new -connections and add them to the network as we did before. +connections and add them to the network as we did before. Do not forget to add +the :code:`x=0` specification to make liquid saturated state at the outlet of +the condenser. The valve connects to the drum at the inlet :code:`'in1'`. The drum's outlet :code:`'out1'` is saturated liquid and connects to the evaporator's cold side @@ -264,16 +266,16 @@ connections to the model: Parametrization +++++++++++++++ Previous parametrization stays untouched. Regarding the evaporator, we specify -pressure ratios on hot side as well as the evaporation pressure, for which we -can obtain a good initial guess based on the ambient temperature level using -CoolProp. From this specification the pinch point layout will be a result, -similar as in waste heat steam generators. The pressure ratio of the cold side -*MUST NOT* be specified in this setup as the drum assumes pressure equality -for all inlets and outlets. +pressure ratios on hot side. On top, the evaporation pressure is indirectly +defined through the evaporation temperature and the fact, that we are in +two-phase region. From this specification the pinch point layout will be a +result, similar as in waste heat steam generators. The pressure ratio of the +cold side *MUST NOT* be specified in this setup as the drum assumes pressure +equality for all inlets and outlets. The superheater will also use the pressure ratios on hot and cold side. -Further we set a value for the enthalpy at the working fluid side outlet. This -determines the degree of overheating and is again based on a good guess. +Further we set a value for the degree of superheating :code:`td_dew` at the +working fluid side outlet. .. literalinclude:: /../tutorial/advanced/stepwise.py :language: python @@ -308,8 +310,8 @@ Compressor system To complete the heat pump, we will add the compressor system to our existing network. This requires to change the connections 0, 6 and 17. The connection 6 has to be changed to include the compressor. After the last compressor stage, -connection 0 has to redefined, since we need to include the CycleCloser of the -working fluid's cycle. The connection 17 has to be connected to the heat +connection 0 has to be redefined, since we need to include the CycleCloser of +the working fluid's cycle. The connection 17 has to be connected to the heat exchanger for intermittent cooling as well as the bypass. .. figure:: /_static/images/tutorials/heat_pump_stepwise/flowsheet.svg @@ -375,18 +377,21 @@ pressure losses on both sides. :start-after: [sec_14] :end-before: [sec_15] -Regarding the connections we set enthalpy values for all working fluid side -connections. After the superheater and intermittent cooling the value will be -near saturation (enthalpy value of connection c5), after the compressors it -will be higher. +Regarding the connections we set degree of superheating at the outlet of the +superheater and intermittent cooling, since these values will be near +saturation. The temperature at both compressor exits will be higher than those, +we can implicitly define this through the enthalpy fixed to the saturated gas +enthalpy multiplied by a bit. This specification is for numerical stability, it +will be changed in the next step. For the ambient side, we set temperature, pressure and fluid at connection 11. On top of that, we can specify the temperature of the ambient water after leaving the intermittent cooler. -With re-adding of connection 0 we have to set the fluid and the pressure again, -but not the temperature value, because this value will be a result of the -condensation pressure and the given enthalpy at the compressor's outlet. +With re-adding of connection 0 we have to set the fluid and the dew line +temperature again, but not the actual temperature value, because this value +will be a result of the condensation pressure and the given enthalpy at the +compressor's outlet. .. literalinclude:: /../tutorial/advanced/stepwise.py :language: python @@ -398,18 +403,16 @@ Solve and Set Final System Parameters Now we solve again. After that, we can exchange our good guesses with actual useful parameters: -The condensation and evaporation pressure levels will be replaced by terminal -temperature values of the condenser and the evaporator respectively. The lower -terminal temperature value of the evaporator :code:`ttd_l` defines the pinch -point. The upper terminal temperature value :code:`ttd_u` of the condenser -defines the condensation pressure. +The condensation dew line temperature and evaporation temperature will be +replaced by pinch temperature difference of the condenser and terminal +temperature difference of the evaporator respectively. The lower terminal +temperature value of the evaporator :code:`ttd_l` is identical to its pinch +point temperature difference. The degree of superheating in the superheater will be determined by the upper -terminal temperature instead of the enthalpy value at connection 6. The outlet +terminal temperature instead of the superheating specification. The outlet enthalpies after both compressors are replaced by the isentropic efficiency -values. Finally, the enthalpy after the intermittent cooling is replaced by -the temperature difference to the boiling point. With this we can ensure, the -working fluid does not start to condensate at the intermittent cooler. +values. .. literalinclude:: /../tutorial/advanced/stepwise.py :language: python @@ -439,8 +442,9 @@ The changes we want to apply can be summarized as follows: - All heat exchangers should be calculated based on their heat transfer coefficient with a characteristic for correction of that value depending - on the change of mass flow (:code:`kA_char`). Therefore, terminal temperature - value specifications need to be added to the design parameters. Also, the + on the change of mass flow (:code:`kA_char` and :code:`UA_cecchinato` for the + condenser). Therefore, pinch temperature difference and terminal temperature + difference value specifications need to be design parameters. Also, the temperature at connection 14 cannot be specified anymore, since it will be a result of the intermittent cooler's characteristics. - Pumps and compressors will have a characteristic function for their @@ -454,7 +458,8 @@ On top of that, for the evaporator the characteristic function of the heat transfer coefficient should follow different data than the default characteristic. The name of that line is 'EVAPORATING FLUID' for the cold side. The default line 'DEFAULT' will be kept for the hot side. These lines -are available in the :ref:`tespy.data ` module. +are available in the :ref:`tespy.data ` module. For the condenser +the specifications are inspired from this paper :cite:`cecchinato2010`. .. attention:: diff --git a/docs/community/community.rst b/docs/community/community.rst index 83bc07f5d..9220c3a0f 100644 --- a/docs/community/community.rst +++ b/docs/community/community.rst @@ -15,17 +15,14 @@ about using the software, you are invited to start a discussion there. Monthly online meeting ====================== -There is a monthly online meeting for all TESPy users and developers. You are -invited to join us on **every 3rd Monday of a month at 17:00 CE(S)T** for a -casual get together. The intent of this meeting is to establish a more active -and well connected network of TESPy users and developers. You can bring your -challenges, questions and ideas and we will discuss them together or plan -potential follow-up meetings. +There are online QA meetings every month, please check the +`oemof calendar `__ for date, time and meeting +link. You are invited to join us with your questions, issues and suggestions! +Please also check the announcement banner at the top of the documentation +pages, you will get updates on rescheduled date/time. -If you are interested, you can simply join the meeting at -https://meet.jit.si/tespy_user_meeting. We are looking forward to seeing you! - -.. attention:: - - Please also check the announcement banner at the top of the documentation - pages, you will get updates on rescheduled date/time. +In-person meetings +================== +Furthermore, there are in-person community meetings. These are held once or +twice a year. To learn about upcoming meetings, follow the blog over at +``__ for the announcements. diff --git a/docs/conf.py b/docs/conf.py index a22cd87b5..497107600 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -435,8 +435,8 @@ def setup(app): } html_theme_options = { - "light_logo": "./images/logo_tespy_mid_christmas.svg", - "dark_logo": "./images/logo_tespy_mid_darkmode_christmas.svg", + "light_logo": "./images/logo_tespy_mid.svg", + "dark_logo": "./images/logo_tespy_mid_darkmode.svg", "announcement": """
""", 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..c33d6926b 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 @@ -52,6 +74,8 @@ Other changes are now named appropriately with a leading underscore. On top, some methods have been renamed to clarify what they do (`PR #875 `__). +- The logging level of some messages has been changed from warning to debug + (`PR #893 `__). Bug Fixes ######### @@ -63,7 +87,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/components/combustion/base.py b/src/tespy/components/combustion/base.py index da4e67a4f..2bc78664e 100644 --- a/src/tespy/components/combustion/base.py +++ b/src/tespy/components/combustion/base.py @@ -993,9 +993,17 @@ def convergence_check(self): if outl.m.val_SI < 0 and outl.m.is_var: outl.m.set_reference_val_SI(10) - if not outl.good_starting_values: - if outl.h.val_SI < 7.5e5 and outl.h.is_var: - outl.h.set_reference_val_SI(1e6) + if not outl.good_starting_values and outl.h.is_var: + if outl.calc_T() < 800: + h = h_mix_pT( + outl.p.val_SI, 800, outl.fluid_data, outl.mixing_rule + ) + outl.h.set_reference_val_SI(h) + elif outl.calc_T() > 1500: + h = h_mix_pT( + outl.p.val_SI, 1500, outl.fluid_data, outl.mixing_rule + ) + outl.h.set_reference_val_SI(h) ###################################################################### # additional checks for performance improvement diff --git a/src/tespy/components/component.py b/src/tespy/components/component.py index 830cd8551..d5a40832a 100644 --- a/src/tespy/components/component.py +++ b/src/tespy/components/component.py @@ -448,12 +448,13 @@ def _setup_user_imposed_constraints(self, row_idx, sum_eq): elif data.is_set: msg = ( - 'All parameters of the component group have to be ' - 'specified! This component group uses the following ' - f'parameters: {", ".join(data.elements)} at ' - f'{self.label}. Group will be set to False.' + f"Not all parameters of the component group {key} " + f"of component {self.label} have to been specified, " + "the group equation will not be applied. This " + "component group uses the following " + f"parameters: {', '.join(data.elements)}." ) - logger.warning(msg) + logger.debug(msg) data.set_attr(is_set=False) else: data.set_attr(is_set=False) diff --git a/src/tespy/components/heat_exchangers/base.py b/src/tespy/components/heat_exchangers/base.py index 54318952d..d565b912b 100644 --- a/src/tespy/components/heat_exchangers/base.py +++ b/src/tespy/components/heat_exchangers/base.py @@ -465,9 +465,9 @@ def calculate_td_log(self): T_o2 = o2.calc_T() if T_i1 <= T_o2: - T_i1 = T_o2 + 0.01 + T_o2 = T_i1 - 0.1 if T_o1 <= T_i2: - T_o1 = T_i2 + 0.01 + T_o1 = T_i2 + 0.1 ttd_u = T_i1 - T_o2 ttd_l = T_o1 - T_i2 @@ -563,8 +563,17 @@ def kA_char_func(self): """ p1 = self.kA_char1.param p2 = self.kA_char2.param - f1 = self.get_char_expr(p1, **self.kA_char1.char_params) - f2 = self.get_char_expr(p2, **self.kA_char2.char_params) + if self.local_offdesign: + design_value = self._connection_offdesign[self.inl[0].label][p1] + actual_value = getattr(self.inl[0], p1).val_SI + f1 = actual_value / design_value + + design_value = self._connection_offdesign[self.inl[1].label][p2] + actual_value = getattr(self.inl[1], p2).val_SI + f2 = actual_value / design_value + else: + f1 = self.get_char_expr(p1, **self.kA_char1.char_params) + f2 = self.get_char_expr(p2, **self.kA_char2.char_params) fkA1 = self.kA_char1.char_func.evaluate(f1) fkA2 = self.kA_char2.char_func.evaluate(f2) @@ -1019,7 +1028,7 @@ def calc_parameters(self): "because cold side inlet temperature is out of bounds for hot " "side fluid." ) - logger.warning(msg) + logger.debug(msg) try: self.eff_cold.val_SI = ( (self.outl[1].h.val_SI - self.inl[1].h.val_SI) @@ -1032,7 +1041,7 @@ def calc_parameters(self): "because hot side inlet temperature is out of bounds for cold " "side fluid." ) - logger.warning(msg) + logger.debug(msg) self.eff_max.val_SI = max(self.eff_hot.val_SI, self.eff_cold.val_SI) def entropy_balance(self): diff --git a/src/tespy/components/heat_exchangers/condenser.py b/src/tespy/components/heat_exchangers/condenser.py index b4cbc1a2a..c39a3ba1b 100644 --- a/src/tespy/components/heat_exchangers/condenser.py +++ b/src/tespy/components/heat_exchangers/condenser.py @@ -421,7 +421,7 @@ def calc_parameters(self): "because cold side inlet temperature is out of bounds for hot " "side fluid." ) - logger.warning(msg) + logger.debug(msg) try: self.eff_cold.val_SI = ( (self.outl[1].h.val_SI - self.inl[1].h.val_SI) @@ -434,7 +434,7 @@ def calc_parameters(self): "because hot side inlet temperature is out of bounds for cold " "side fluid." ) - logger.warning(msg) + logger.debug(msg) self.eff_max.val_SI = max(self.eff_hot.val_SI, self.eff_cold.val_SI) def convergence_check(self): diff --git a/src/tespy/components/heat_exchangers/sectioned.py b/src/tespy/components/heat_exchangers/sectioned.py index 082961a9f..58bdfe7b2 100644 --- a/src/tespy/components/heat_exchangers/sectioned.py +++ b/src/tespy/components/heat_exchangers/sectioned.py @@ -759,9 +759,9 @@ def UA_cecchinato_func(self): secondary_index = 0 m_r = self.inl[refrigerant_index].m - m_ratio_r = m_r.val_SI / m_r.design + m_ratio_r = max(m_r.val_SI / m_r.design, 1e-6) m_sf = self.inl[secondary_index].m - m_ratio_sf = m_sf.val_SI / m_sf.design + m_ratio_sf = max(m_sf.val_SI / m_sf.design, 1e-6) fUA = ( (1 + alpha_ratio * area_ratio) diff --git a/src/tespy/components/nodes/humidity_control.py b/src/tespy/components/nodes/humidity_control.py new file mode 100644 index 000000000..1d63d4279 --- /dev/null +++ b/src/tespy/components/nodes/humidity_control.py @@ -0,0 +1,324 @@ +# -*- coding: utf-8 + +"""Module of class HumidityControl. + + +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/components/nodes/humidity_control.py + +SPDX-License-Identifier: MIT +""" + +from tespy.components.component import component_registry +from tespy.components.nodes.base import NodeBase +from tespy.tools.data_containers import ComponentMandatoryConstraints as dc_cmc +from tespy.tools.data_containers import SimpleDataContainer as dc_simple +from tespy.tools.fluid_properties import dT_mix_dph +from tespy.tools.fluid_properties import dT_mix_pdh + +# from tespy.tools.fluid_properties import dT_mix_ph_dfluid + + +@component_registry +class HumidityControl(NodeBase): + r""" + A humidity control handles the water content from a humid air flow. + + **Mandatory Equations** + + - :py:meth:`tespy.components.nodes.base.NodeBase.mass_flow_func` + - :py:meth:`tespy.components.nodes.base.NodeBase.pressure_structure_matrix` + - :py:meth:`tespy.components.nodes.humidity_control.HumidityControl.fluid_func` + - :py:meth:`tespy.components.nodes.humidity_control.HumidityControl.energy_balance_func` + + Inlets/Outlets + + - in1: inlet with humid air + - out1: outlet with humid air + - out2: outlet of liquid water or ice + + Image + + .. image:: /api/_images/Splitter.svg + :alt: flowsheet of the splitter + :align: center + :class: only-light + + .. image:: /api/_images/Splitter_darkmode.svg + :alt: flowsheet of the splitter + :align: center + :class: only-dark + + Note + ---- + Fluid separation requires power and cooling, equations have not been + implemented, yet! + + Parameters + ---------- + label : str + The label of the component. + + design : list + List containing design parameters (stated as String). + + offdesign : list + List containing offdesign parameters (stated as String). + + design_path : str + Path to the components design case. + + local_offdesign : boolean + Treat this component in offdesign mode in a design calculation. + + local_design : boolean + Treat this component in design mode in an offdesign calculation. + + char_warnings : boolean + Ignore warnings on default characteristics usage for this component. + + printout : boolean + Include this component in the network's results printout. + + num_out : float, dict + Number of outlets for this component, default value: 2. + + Example + ------- + The separator is used to split up a single mass flow into a specified + number of different parts at identical pressure and temperature but + different fluid composition. Fluids can be separated from each other. + + >>> from tespy.components import Sink, Source, Separator + >>> from tespy.connections import Connection + >>> from tespy.networks import Network + >>> nw = Network(iterinfo=False) + >>> nw.units.set_defaults(**{ + ... "pressure": "bar", "temperature": "degC" + ... }) + >>> so = Source('source') + >>> si1 = Sink('sink1') + >>> si2 = Sink('sink2') + >>> s = Separator('separator', num_out=2) + >>> inc = Connection(so, 'out1', s, 'in1') + >>> outg1 = Connection(s, 'out1', si1, 'in1') + >>> outg2 = Connection(s, 'out2', si2, 'in1') + >>> nw.add_conns(inc, outg1, outg2) + + An Air (simplified) mass flow of 5 kg/s is split up into two mass flows. + One mass flow of 1 kg/s containing 10 % oxygen and 90 % nitrogen leaves the + separator. It is possible to calculate the fluid composition of the second + mass flow. Specify starting values for the second mass flow fluid + composition for calculation stability. + + >>> inc.set_attr(fluid={'O2': 0.23, 'N2': 0.77}, p=1, T=20, m=5) + >>> outg1.set_attr(fluid={'O2': 0.1, 'N2': 0.9}, m=1) + >>> outg2.set_attr(fluid0={'O2': 0.5, 'N2': 0.5}) + >>> nw.solve('design') + >>> outg2.fluid.val['O2'] + 0.2625 + + In the same way, it is possible to specify one of the fluid components in + the second mass flow instead of the first mass flow. The solver will find + the mass flows matching the desired composition. 65 % of the mass flow + will leave the separator at the second outlet the case of 30 % oxygen + mass fraction for this outlet. + + >>> outg1.set_attr(m=None) + >>> outg2.set_attr(fluid={'O2': 0.3}) + >>> nw.solve('design') + >>> outg2.fluid.val['O2'] + 0.3 + >>> round(outg2.m.val_SI / inc.m.val_SI, 2) + 0.65 + """ + + @staticmethod + def get_parameters(): + return {'num_out': dc_simple(description="number of outlets")} + + def _update_num_eq(self): + self.variable_fluids = set( + [fluid for c in self.inl + self.outl for fluid in c.fluid.is_var] + ) + num_fluid_eq = len(self.variable_fluids) + if num_fluid_eq == 0: + num_fluid_eq = 1 + self.variable_fluids = [list(self.inl[0].fluid.is_set)[0]] + + self.constraints["fluid_constraints"].num_eq = num_fluid_eq + + def get_mandatory_constraints(self): + return { + 'mass_flow_constraints': dc_cmc(**{ + 'num_eq_sets': 1, + 'func': self.mass_flow_func, + 'dependents': self.mass_flow_dependents, + 'description': 'mass balance constraint' + }), + 'fluid_constraints': dc_cmc(**{ + 'num_eq_sets': self.num_o, + 'func': self.fluid_func, + 'deriv': self.fluid_deriv, + 'dependents': self.fluid_dependents, + 'description': 'fluid mass fraction balance constraints' + }), + 'energy_balance_constraints': dc_cmc(**{ + 'num_eq_sets': self.num_o, + 'func': self.energy_balance_func, + 'deriv': self.energy_balance_deriv, + 'dependents': self.energy_balance_dependents, + 'description': 'equal temperature at all outlets constraints' + }), + 'pressure_constraints': dc_cmc(**{ + 'num_eq_sets': self.num_o, + 'structure_matrix': self.pressure_structure_matrix, + 'description': 'pressure equality constraints' + }) + } + + @staticmethod + def inlets(): + return ['in1'] + + def outlets(self): + if self.num_out.is_set: + return [f'out{i + 1}' for i in range(self.num_out.val)] + else: + self.set_attr(num_out=2) + return self.outlets() + + def propagate_wrapper_to_target(self, branch): + branch["components"] += [self] + for outconn in self.outl: + branch["connections"] += [outconn] + outconn.target.propagate_wrapper_to_target(branch) + + def fluid_func(self): + r""" + Calculate the vector of residual values for fluid balance equations. + + Returns + ------- + residual : list + Vector of residual values for component's fluid balance. + + .. math:: + + 0 = \dot{m}_{in} \cdot x_{fl,in} - \dot {m}_{out,j} + \cdot x_{fl,out,j}\\ + \forall fl \in \text{network fluids,} + \; \forall j \in \text{outlets} + """ + i = self.inl[0] + residual = [] + for fluid in self.variable_fluids: + res = i.fluid.val[fluid] * i.m.val_SI + for o in self.outl: + res -= o.fluid.val[fluid] * o.m.val_SI + residual += [res] + return residual + + def fluid_deriv(self, increment_filter, k, dependents=None): + r""" + Calculate partial derivatives of fluid balance. + + Parameters + ---------- + increment_filter : ndarray + Matrix for filtering non-changing variables. + + k : int + Position of derivatives in Jacobian matrix (k-th equation). + """ + i = self.inl[0] + for fluid in self.variable_fluids: + for o in self.outl: + self._partial_derivative(o.m, k, -o.fluid.val[fluid], increment_filter) + if fluid in o.fluid.is_var: + self.jacobian[k, o.fluid.J_col[fluid]] = -o.m.val_SI + + self._partial_derivative(i.m, k, i.fluid.val[fluid], increment_filter) + if fluid in i.fluid.is_var: + self.jacobian[k, i.fluid.J_col[fluid]] = i.m.val_SI + + k += 1 + + def fluid_dependents(self): + return { + "scalars": [ + [c.m for c in self.inl + self.outl] + for f in self.variable_fluids + ], + "vectors": [{ + c.fluid: set(f) & c.fluid.is_var for c in self.inl + self.outl + } for f in self.variable_fluids] + } + + def energy_balance_func(self): + r""" + Calculate energy balance. + + Returns + ------- + residual : list + Residual value of energy balance. + + .. math:: + + 0 = T_{in} - T_{out,j}\\ + \forall j \in \text{outlets} + """ + residual = [] + T_in = self.inl[0].calc_T() + for o in self.outl: + residual += [T_in - o.calc_T()] + return residual + + def energy_balance_deriv(self, increment_filter, k, dependents=None): + r""" + Calculate partial derivatives of energy balance. + + Parameters + ---------- + increment_filter : ndarray + Matrix for filtering non-changing variables. + + k : int + Position of derivatives in Jacobian matrix (k-th equation). + """ + i = self.inl[0] + dT_dp_in = 0 + dT_dh_in = 0 + if i.p.is_var: + # outlet pressure must be variable as well in this case! + dT_dp_in = dT_mix_dph(i.p.val_SI, i.h.val_SI, i.fluid_data, i.mixing_rule) + if i.h.is_var: + dT_dh_in = dT_mix_pdh(i.p.val_SI, i.h.val_SI, i.fluid_data, i.mixing_rule) + + for o in self.outl: + args = (o.p.val_SI, o.h.val_SI, o.fluid_data, o.mixing_rule) + + dT_dp_out = 0 + if o.p.is_var: + dT_dp_out = -dT_mix_dph(*args) + # pressure is always coupled + self._partial_derivative(i.p, k, dT_dp_in - dT_dp_out) + + if o.h.is_var: + dT_dh_out = -dT_mix_pdh(*args) + + # enthalpy is not necessarily coupled + if i.h._reference_container == o.h._reference_container: + self._partial_derivative(i.h, k, dT_dh_in - dT_dh_out) + else: + self._partial_derivative(i.h, k, dT_dh_in) + self._partial_derivative(o.h, k, dT_dh_out) + + k += 1 + + def energy_balance_dependents(self): + return [ + [self.inl[0].p, self.inl[0].h, o.p, o.h] for o in self.outl + ] diff --git a/src/tespy/connections/__init__.py b/src/tespy/connections/__init__.py index 3741bb556..2eba79104 100644 --- a/src/tespy/connections/__init__.py +++ b/src/tespy/connections/__init__.py @@ -3,4 +3,5 @@ from .bus import Bus # noqa: F401 from .connection import Connection # noqa: F401 from .connection import Ref # noqa: F401 +from .humidairconnection import HAConnection # noqa: F401 from .powerconnection import PowerConnection # noqa: F401 diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index e48962192..76eee5f24 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -56,6 +56,7 @@ from tespy.tools.helpers import _is_variable from tespy.tools.helpers import _partial_derivative from tespy.tools.helpers import _partial_derivative_vecvar +from tespy.tools.helpers import seeded_random from tespy.tools.units import SI_UNITS @@ -147,6 +148,23 @@ def _parameter_specification(self, key, value): elif is_numeric or is_quantity: # value specification if key in self.property_data: + if f"{key}_ref" in self.property_data: + if self.get_attr(f"{key}_ref" ).is_set: + msg = ( + f"You have specified a numerical value for the " + f"parameter '{key}' while having a Ref specified " + f"at the same time at connection {self.label}. In " + "the moment, this does not overwrite setting the " + "Ref. To unset the Ref before setting the " + f"numerical value, run .set_attr({key}=None) " + f"before .set_attr({key}={value}). With the next " + "major release of TESPy setting a numerical value " + "will always replace a previously specified Ref " + "making it impossible to set a numerical value " + "and a Ref for one parameter on one connection " + "simultaneously." + ) + warnings.warn(msg, FutureWarning) self.get_attr(key).set_attr(is_set=True, val=value) else: self.get_attr(key.replace('0', '')).set_attr(val0=value) @@ -158,6 +176,22 @@ def _parameter_specification(self, key, value): logger.error(msg) raise NotImplementedError(msg) else: + if self.get_attr(key).is_set: + msg = ( + f"You have specified a Ref for the parameter '{key}' " + "while having a numerical value specified at the same " + f"time at connection {self.label}. In the moment, " + "this does not overwrite setting the numerical value. " + "To unset the specified value before setting the Ref, " + f"run .set_attr({key}=None) before " + f".set_attr({key}=Ref(...)). With the next major " + "release of TESPy setting a Ref will automatically " + "replace a previously specified numerical value " + "making it impossible to set a numerical value and a " + "Ref for one parameter on one connection " + "simultaneously." + ) + warnings.warn(msg, FutureWarning) self.get_attr(f"{key}_ref").set_attr(ref=value) self.get_attr(f"{key}_ref").set_attr(is_set=True) @@ -212,6 +246,9 @@ def _serializable(): def get_variables(self): return {} + def _guess_starting_values(self, units): + pass + def _preprocess(self, row_idx): self.num_eq = 0 @@ -724,7 +761,7 @@ def set_attr(self, **kwargs): # other boolean keywords elif key in ['printout', 'local_design', 'local_offdesign']: if not isinstance(kwargs[key], bool): - msg = ('Please provide the ' + key + ' as boolean.') + msg = f"Please provide the {key} as boolean." logger.error(msg) raise TypeError(msg) else: @@ -772,6 +809,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 +873,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 +883,79 @@ 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 _guess_starting_values(self, units): + # the below part does not work for PowerConnection right now + if sum(self.fluid.val.values()) == 0: + msg = ( + 'The starting value for the fluid composition of the ' + f'connection {self.label} is empty. This might lead to issues ' + 'in the initialisation and solving process as fluid ' + 'property functions can not be called. Make sure you ' + 'specified a fluid composition in all parts of the network.' + ) + logger.warning(msg) + + for key, variable in self.get_variables().items(): + # for connections variables can be presolved and not be var anymore + if variable.is_var: + if not self.good_starting_values: + self._guess_starting_value_from_connected_components(key, units) + + variable.set_SI_from_val0(units) + # variable.set_SI_from_val0() + variable.set_reference_val_SI(variable._val_SI) + + self._precalc_guess_values() + + def _guess_starting_value_from_connected_components(self, key, units): + r""" + Set starting values for fluid properties. + + The component classes provide generic starting values for their inlets + and outlets. + + Parameters + ---------- + c : tespy.connections.connection.Connection + Connection to initialise. + """ + if np.isnan(self.get_attr(key).val0): + # starting value for mass flow is random between 1 and 2 kg/s + # (should be generated based on some hash maybe?) + if key == 'm': + rndm = seeded_random(self.label) + value = float(rndm + 1) + + # generic starting values for pressure and enthalpy + elif key in ['p', 'h']: + # retrieve starting values from component information + val_s = self.source.initialise_source(self, key) + val_t = self.target.initialise_target(self, key) + + if val_s == 0 and val_t == 0: + if key == 'p': + value = 1e5 + elif key == 'h': + value = 1e6 + + elif val_s == 0: + value = val_t + elif val_t == 0: + value = val_s + else: + value = (val_s + val_t) / 2 + + # these values are SI, so they are set to the respective variable + self.get_attr(key).set_reference_val_SI(value) + self.get_attr(key).set_val0_from_SI(units) def _precalc_guess_values(self): """ @@ -926,6 +1039,7 @@ def _presolve(self): raise TESPyNetworkError(msg) presolved_equations = [] + if self.p.is_set: if self.T_dew.is_set or self.T_bubble.is_set: msg = ( @@ -1197,7 +1311,7 @@ def get_parameters(self): dependents=self.Td_bp_dependents, num_eq=1, quantity="temperature_difference", description="temperature difference to boiling point (deprecated)" - ) + ), } def get_fluid_data(self): @@ -1719,7 +1833,7 @@ def _adjust_to_property_limits(self, nw): self._adjust_enthalpy(fl) # two-phase related - if (self.Td_bp.is_set or self.state.is_set or self.x.is_set or self.td_bubble.is_set or self.td_dew.is_set) and self.it < 10: + if (self.Td_bp.is_set or self.state.is_set or self.x.is_set or self.td_bubble.is_set or self.td_dew.is_set) and self.it < 30: self._adjust_to_two_phase(fl) # mixture diff --git a/src/tespy/connections/humidairconnection.py b/src/tespy/connections/humidairconnection.py new file mode 100644 index 000000000..807196845 --- /dev/null +++ b/src/tespy/connections/humidairconnection.py @@ -0,0 +1,260 @@ +# -*- coding: utf-8 +"""Module of class Connection and class Ref. + + +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/connections/humidairconnection.py + +SPDX-License-Identifier: MIT +""" + +import numpy as np +from CoolProp.CoolProp import HAPropsSI + +from tespy.tools import fluid_properties as fp +from tespy.tools.data_containers import FluidComposition as dc_flu +from tespy.tools.data_containers import FluidProperties as dc_prop +from tespy.tools.data_containers import SimpleDataContainer as dc_simple +from tespy.tools.fluid_properties.functions import w_mix_pT_humidair +from tespy.tools.fluid_properties.functions import h_mix_pT +from tespy.tools.fluid_properties.mixtures import _get_fluid_alias, w_mix_fluid_data +from tespy.tools.helpers import seeded_random + +from .connection import Connection +from .connection import connection_registry + + +@connection_registry +class HAConnection(Connection): + + def get_parameters(self): + return { + "fluid": dc_flu( + d=1e-5, + description="mass fractions of the fluid composition" + ), + "m": dc_prop( + quantity="mass_flow", + description="mass flow of dry air (system variable)" + ), + "mHA": dc_prop( + quantity="mass_flow", + description="mass flow of humid air" + ), + "mH2O": dc_prop( + quantity="mass_flow", + description="mass flow of liquid or solid water not contained in humid air" + ), + "p": dc_prop( + quantity="pressure", + description="absolute pressure of the fluid (system variable)" + ), + "h": dc_prop( + quantity="enthalpy", + description="dry air mass specific enthalpy (system variable)" + ), + "w": dc_prop( + quantity="ratio", + description="mass of water per mass of dry air" + ), + "T": dc_prop( + func=self.T_func, + dependents=self.T_dependents, + num_eq=1, + quantity="temperature", + description="temperature of the fluid" + ), + "v": dc_prop( + func=self.v_func, + dependents=self.v_dependents, + num_eq=1, + quantity="volumetric_flow", + description="volumetric flow of the fluid" + ), + "vol": dc_prop( + quantity="specific_volume", + description="specific volume of the fluid (output only)" + ), + "s": dc_prop( + quantity="entropy", + description="specific entropy of the fluid (output only)" + ), + "r": dc_prop( + func=self.r_func, + dependents=self.r_dependents, + num_eq=1, + quantity="ratio", + description="relative humidity" + ), + "fluid_balance": dc_simple( + func=self.fluid_balance_func, + deriv=self.fluid_balance_deriv, + _val=False, num_eq_sets=1, + dependents=self.fluid_balance_dependents, + description="apply an equation which closes the fluid balance with at least two unknown fluid mass fractions" + ), + } + + def _parameter_specification(self, key, value): + if key == "w" or key == "w0": + # specification of w is equivalent to specification of fluid + # composition for humid air + air = 1 / (1 + value) + if key == "w": + self.set_attr(fluid={"air": air, "water": 1 - air}) + else: + self.set_attr(fluid0={"air": air, "water": 1 - air}) + else: + super()._parameter_specification(key, value) + + # for HAConnection mixing rule cannot be modified, is always humidair + def _get_mixing_rule(self): + return "humidair" + + def _set_mixing_rule(self, value): + if value is not None and value != self.mixing_rule: + print(value) + msg = ( + "You cannot change the mixing rule specification for a " + f"Connection of type {self.__class__.__name__}" + ) + raise ValueError(msg) + + mixing_rule = property(_get_mixing_rule, _set_mixing_rule) + + def _guess_starting_values(self, units): + if self.h.is_var and not self.good_starting_values: + value = seeded_random(self.label) + T_rand = 280 + value * (300 - 280) + h = fp.h_mix_pT(1e5, T_rand, self.fluid_data, self.mixing_rule) + self.h.set_reference_val_SI(h) + self._precalc_guess_values() + + def _precalc_guess_values(self): + if not self.h.is_var: + return + + if not self.good_starting_values: + if self.T.is_set: + try: + w = self.calc_w() + self.h.set_reference_val_SI( + HAPropsSI("H", "P", self.p.val_SI, "T", self.T.val_SI, "W", w) + ) + except ValueError: + pass + + def _presolve(self): + + air_alias = _get_fluid_alias("air", self.fluid_data) + water_alias = _get_fluid_alias("water", self.fluid_data) + if not air_alias: + msg = "air must be present in fluid composition" + raise ValueError(msg) + + elif not water_alias: + msg = "water must be present in fluid composition" + raise ValueError(msg) + + if len(self.fluid.is_var) > 0: + return [] + + presolved_equations = [] + + if self.h.is_var and not self.p.is_var: + if self.T.is_set: + self.h.set_reference_val_SI(h_mix_pT(self.p.val_SI, self.T.val_SI, self.fluid_data, self.mixing_rule)) + self.h._potential_var = False + if "T" in self._equation_set_lookup.values(): + presolved_equations += ["T"] + + presolved_equations = [ + key for parameter in presolved_equations + for key, value in self._equation_set_lookup.items() + if value == parameter + ] + return presolved_equations + + def _adjust_to_property_limits(self, nw): + + if self.p.is_var: + if self.p.val_SI < 100: + self.p.val_SI = 101 + elif self.p.val_SI > 100e5: + self.p.val_SI = 99e5 + + if self.h.is_var: + # TODO: check minimum temperature how it matches minimum humidity ratio + d = self.h._reference_container._d + hmin = HAPropsSI("H", "T", -50 + 273.15, "P", self.p.val_SI, "R", 1) + if self.h.val_SI < hmin: + delta = max(abs(self.h.val_SI * d), d) * 5 + self.h.set_reference_val_SI(hmin + delta) + + else: + # TODO: where to get reasonable hmax from?! + hmax = HAPropsSI("H", "T", 300 + 273.15, "P", self.p.val_SI, "R", 0) + if self.h.val_SI > hmax: + delta = max(abs(self.h.val_SI * d), d) * 5 + self.h.set_reference_val_SI(hmax - delta) + + @classmethod + def _result_attributes(cls): + return ["m", "mHA", "mH2O", "p", "h", "T", "w", "s", "vol", "v", "r"] + + @classmethod + def _print_attributes(cls): + return ["m", "mHA", "mH2O", "p", "h", "T", "w", "r"] + + def calc_r(self): + w = self.calc_w() + try: + return HAPropsSI("R", "P", self.p.val_SI, "T", self.T.val_SI, "W", w) + except ValueError as e: + value = str(e).split("value (")[1].split(")")[0] + return float(value) + + def r_func(self): + return self.r.val_SI - self.calc_r() + + def r_dependents(self): + water_alias = _get_fluid_alias("H2O", self.fluid_data) + # water alias is already a set + return { + "scalars": [self.p, self.h], + "vectors": [{self.fluid: water_alias}] + } + + def calc_w(self): + return w_mix_pT_humidair(self.p.val_SI, self.T.val_SI, self.fluid_data) + + def calc_results(self, units): + self.T.val_SI = self.calc_T() + self.vol.val_SI = self.calc_vol() # Mixture volume per mass of dry air + self.v.val_SI = self.vol.val_SI * self.m.val_SI # Mixture volume flow rate + # handle the water fraction + self.w.val_SI = self.calc_w() + # # Convert from kg water/kg dry air to mass fraction of water in humid air + # x_h2o = self.w.val_SI / (1 + self.w.val_SI) + # x_air = 1 - x_h2o + self.mHA.val_SI = self.m.val_SI * (1 + self.w.val_SI) + w_mixture = w_mix_fluid_data(self.fluid_data) + # Calculate kg water/kg dry air that is not in the humid air + delta_w = w_mixture - self.w.val_SI + self.mH2O.val_SI = self.m.val_SI * delta_w + self.r.val_SI = self.calc_r() + # if self.r.val_SI > 1: + # self.r.val_SI = np.nan + + for prop in self._result_attributes(): + param = self.get_attr(prop) + if not param.is_set: + param.set_val_from_SI(units) + + self.m.set_val0_from_SI(units) + self.p.set_val0_from_SI(units) + self.h.set_val0_from_SI(units) + self.fluid.val0 = self.fluid.val.copy() + + return True diff --git a/src/tespy/connections/powerconnection.py b/src/tespy/connections/powerconnection.py index 3503e346a..7589c0a41 100644 --- a/src/tespy/connections/powerconnection.py +++ b/src/tespy/connections/powerconnection.py @@ -130,6 +130,10 @@ def set_attr(self, **kwargs): logger.error(msg) raise KeyError(msg) + def _guess_starting_values(self, units): + if self.E.is_var and not self.good_starting_values: + self.E.set_reference_val_SI(0.0) + def _precalc_guess_values(self): pass diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index ed4ae19d5..d2e574adc 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -16,6 +16,7 @@ import json import math import os +from pathlib import Path import warnings from time import time @@ -917,16 +918,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 +963,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 +987,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() @@ -1757,6 +1765,26 @@ def _prepare_design(self): data = _local_designs[path][c] # write data self._write_design_state_to_component(cp, data) + # this is a hack to write the connection specifications of the + # respective component into the component. In the component + # itself, the value from here should be utilized instead of the + # .design value of the connection in case local_offdesign is + # set to True + if cp.inl[0].label in _local_designs[path][cp.inl[0].__class__.__name__].index: + cp._connection_offdesign = { + c.label: _local_designs[path][c.__class__.__name__].loc[c.label] + for c in cp.inl + cp.outl + } + else: + # remap to actual inlet and outlet labels + cp._connection_offdesign = { + c.label: _local_designs[path][c.__class__.__name__].loc[f"hx_in{k + 1}"] + for k, c in enumerate(cp.inl) + } + cp._connection_offdesign.update({ + c.label: _local_designs[path][c.__class__.__name__].loc[f"hx_out{k + 1}"] + for k, c in enumerate(cp.outl) + }) # unset design parameters for var in cp.design: @@ -2023,27 +2051,7 @@ def _set_starting_values(self): df = dfs[c.__class__.__name__] self._write_starting_values_to_connection(c, df) - if type(c) == Connection: - # the below part does not work for PowerConnection right now - if sum(c.fluid.val.values()) == 0: - msg = ( - 'The starting value for the fluid composition of the ' - f'connection {c.label} is empty. This might lead to issues ' - 'in the initialisation and solving process as fluid ' - 'property functions can not be called. Make sure you ' - 'specified a fluid composition in all parts of the network.' - ) - logger.warning(msg) - - for key, variable in c.get_variables().items(): - # for connections variables can be presolved and not be var anymore - if variable.is_var: - if not c.good_starting_values: - self._guess_starting_value_from_connected_components(c, key) - - variable.set_SI_from_val0(self.units) - # variable.set_SI_from_val0() - variable.set_reference_val_SI(variable._val_SI) + c._guess_starting_values(self.units) for cp in self.comps["object"]: for key, variable in cp.get_variables().items(): @@ -2054,70 +2062,35 @@ def _set_starting_values(self): variable.set_SI_from_val(self.units) variable.set_reference_val_SI(variable._val_SI) - for c in self.conns['object']: - c._precalc_guess_values() - msg = 'Generic fluid property specification complete.' logger.debug(msg) - def _guess_starting_value_from_connected_components(self, c, key): - r""" - Set starting values for fluid properties. - - The component classes provide generic starting values for their inlets - and outlets. - - Parameters - ---------- - c : tespy.connections.connection.Connection - Connection to initialise. - """ - if np.isnan(c.get_attr(key).val0): - # starting value for mass flow is random between 1 and 2 kg/s - # (should be generated based on some hash maybe?) - if key == 'm': - seed = abs(hash(c.label)) % (2**32) - rng = np.random.default_rng(seed=seed) - value = float(rng.random() + 1) - - # generic starting values for pressure and enthalpy - elif key in ['p', 'h']: - # retrieve starting values from component information - val_s = c.source.initialise_source(c, key) - val_t = c.target.initialise_target(c, key) - - if val_s == 0 and val_t == 0: - if key == 'p': - value = 1e5 - elif key == 'h': - value = 1e6 - - elif val_s == 0: - value = val_t - elif val_t == 0: - value = val_s - else: - value = (val_s + val_t) / 2 - - elif key == 'E': - value = 0.0 - - # these values are SI, so they are set to the respective variable - c.get_attr(key).set_reference_val_SI(value) - c.get_attr(key).set_val0_from_SI(self.units) @staticmethod - def _load_network_state(json_path): + def _load_network_state(json_path: str | bytes | bytearray | Path): r""" Read network state from given file. Parameters ---------- - json_path : str + json_path : str | bytes | bytearray | Path Path to network information. """ - with open(json_path, "r") as f: - data = json.load(f) + data = None + if not isinstance(json_path, Path): + try: + data = json.loads(json_path) + except json.JSONDecodeError as e: + msg = ( + "The provided json_path could not be decoded. If this is not " + "a valid json string, please provide a valid file path instead of " + "%s" + ) + logger.debug(msg, str(json_path)) + pass + if data is None: + with open(json_path, "r") as f: + data = json.load(f) dfs = {} if "Connection" in data["Connection"]: @@ -2358,6 +2331,23 @@ def _get_linear_dependents_by_variable_index(self, idx) -> list: variables = [self._variable_lookup[v] for v in dependents["variables"]] variable_list = [(v["object"].label, v["property"]) for v in variables] return variable_list + + def get_sorted_residual_index(self) -> list[int]: + """Get the sorted array of residual indices. + + Returns + ------- + list[int] + List of variable numbers, the index values. + """ + # vars: dict[tuple[int, str], dict] = self.get_variables() + sidx: list[int] = list(np.argsort(np.abs(self.residual))[::-1]) + # sres = np.array([self.residual[i] for i in sidx]) + # chis = self.residual_history.shape[1] + # for i in range(2, n): + # sres = np.vstack((sres, [self.residual_history[i-2][j] for j in sidx])) + # sres = np.vstack((sres, self.residual_history[-n+1:, :][:, sidx].T)) + return sidx def solve(self, mode, init_path=None, design_path=None, max_iter=50, min_iter=4, init_only=False, init_previous=True, @@ -2845,18 +2835,19 @@ def _update_variables(self): # get_J_col yet relax = 1 if self.robust_relax: - if self.iter < 3: - relax = 0.25 - elif self.iter < 5: - relax = 0.5 - elif self.iter < 8: - relax = 0.75 + # relax_values = [ + # (0.05, 0.05 * self.max_iter), + # (0.10, 0.10 * self.max_iter), + # (0.25, 0.25 * self.max_iter), + # (0.50, 0.50 * self.max_iter) + # ] + relax = 0.05 + 0.95 * min(1, self.iter / (0.25 * self.max_iter)) for _, data in self.variables_dict.items(): if data["variable"] in ["m", "h", "E"]: container = data["obj"] container._val_SI += increment[container.J_col] * relax - elif data["variable"] == "p": + elif data["variable"] in ["p", "w"]: container = data["obj"] p_relax = max( 1, -2 * increment[container.J_col] / container.val_SI @@ -2908,7 +2899,7 @@ def _adapt_to_variable_bounds(self): # - only in design case if ( self.iter < 3 - and norm(self.increment) > 1e3 + and norm(self.increment) > 1e-1 and self.mode == "design" ): for cp in self.comps['object']: @@ -3549,18 +3540,26 @@ def export(self, json_file_path=None): return export - def save(self, json_file_path): + def save(self, json_file_path: str | Path | None) -> None | str: r""" Dump the results to a json style output. Parameters ---------- - json_file_path : str + json_file_path : str | Path | None Filename to dump results into. + Returns + ------- + None + If a file path is provided, results are saved to file. + str + If no file path is provided, results are returned as string. + Note ---- - Results will be saved to specified file path + Results will be saved to specified file path in json format. If no + file path is provided, the results will be returned as string. """ dump = {} @@ -3571,6 +3570,9 @@ def save(self, json_file_path): dump = hlp._nested_dict_of_dataframes_to_dict(dump) + if json_file_path is None: + return json.dumps(dump, indent=2) + with open(json_file_path, "w") as f: json.dump(dump, f) 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..7b19da5bd 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/src/tespy/tools/fluid_properties/functions.py @@ -10,6 +10,7 @@ SPDX-License-Identifier: MIT """ +from CoolProp.CoolProp import HAPropsSI from tespy.tools.global_vars import FLUID_ALIASES from tespy.tools.logger import logger @@ -25,6 +26,7 @@ from .mixtures import T_MIX_PS_REVERSE from .mixtures import V_MIX_PT_DIRECT from .mixtures import VISCOSITY_MIX_PT_DIRECT +from .mixtures import w_mix_pT_humidair, w_mix_ph_humidair, w_mix_ps_humidair def isentropic(p_1, h_1, p_2, fluid_data, mixing_rule=None, T0=None): @@ -138,11 +140,15 @@ def T_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): return pure_fluid["wrapper"].T_ph(p, h) else: _check_mixing_rule(mixing_rule, T_MIX_PH_REVERSE, "temperature (from enthalpy)") - kwargs = { - "p": p, "target_value": h, "fluid_data": fluid_data, "T0": T0, - "f": T_MIX_PH_REVERSE[mixing_rule] - } - return inverse_temperature_mixture(**kwargs) + if mixing_rule == "humidair": + w = w_mix_ph_humidair(p, h, fluid_data) + return HAPropsSI("T", "P", p, "H", h, "W", w) + else: + kwargs = { + "p": p, "target_value": h, "fluid_data": fluid_data, "T0": T0, + "f": T_MIX_PH_REVERSE[mixing_rule] + } + return inverse_temperature_mixture(**kwargs) def dT_mix_pdh(p, h, fluid_data, mixing_rule=None, T0=None): @@ -296,12 +302,16 @@ def T_mix_ps(p, s, fluid_data, mixing_rule=None, T0=None): pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].T_ps(p, s) else: - _check_mixing_rule(mixing_rule, T_MIX_PS_REVERSE, "temperature (from entropy)") - kwargs = { - "p": p, "target_value": s, "fluid_data": fluid_data, "T0": T0, - "f": T_MIX_PS_REVERSE[mixing_rule] - } - return inverse_temperature_mixture(**kwargs) + if mixing_rule == "humidair": + w = w_mix_ps_humidair(p, s, fluid_data) + return HAPropsSI("T", "P", p, "S", s, "W", w) + else: + _check_mixing_rule(mixing_rule, T_MIX_PS_REVERSE, "temperature (from entropy)") + kwargs = { + "p": p, "target_value": s, "fluid_data": fluid_data, "T0": T0, + "f": T_MIX_PS_REVERSE[mixing_rule] + } + return inverse_temperature_mixture(**kwargs) def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): @@ -309,8 +319,12 @@ def v_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): pure_fluid = get_pure_fluid(fluid_data) return 1 / pure_fluid["wrapper"].d_ph(p, h) else: - T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) - return v_mix_pT(p, T, fluid_data, mixing_rule) + if mixing_rule == "humidair": + w = w_mix_ph_humidair(p, h, fluid_data) + return HAPropsSI("V", "P", p, "H", h, "W", w) + else: + T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) + return v_mix_pT(p, T, fluid_data, mixing_rule) def dv_mix_dph(p, h, fluid_data, mixing_rule=None, T0=None): @@ -341,8 +355,12 @@ def viscosity_mix_ph(p, h, fluid_data, mixing_rule=None, T0=None): pure_fluid = get_pure_fluid(fluid_data) return pure_fluid["wrapper"].viscosity_ph(p, h) else: - T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) - return viscosity_mix_pT(p, T, fluid_data, mixing_rule) + if mixing_rule == "humidair": + w = w_mix_ph_humidair(p, h, fluid_data) + return HAPropsSI("Visc", "P", p, "H", h, "W", w) + else: + T = T_mix_ph(p, h , fluid_data, mixing_rule, T0) + return viscosity_mix_pT(p, T, fluid_data, mixing_rule) def viscosity_mix_pT(p, T, fluid_data, mixing_rule=None): @@ -352,3 +370,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..7b2f7f8f1 100644 --- a/src/tespy/tools/fluid_properties/helpers.py +++ b/src/tespy/tools/fluid_properties/helpers.py @@ -49,6 +49,9 @@ def get_pure_fluid(fluid_data): def single_fluid(fluid_data): r"""Return the name of the pure fluid in a fluid vector.""" + if "_HUMID_AIR" in fluid_data: + return None + if get_number_of_fluids(fluid_data) > 1: return None else: @@ -360,3 +363,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/mixtures.py b/src/tespy/tools/fluid_properties/mixtures.py index d26cdcb10..8b1a42139 100644 --- a/src/tespy/tools/fluid_properties/mixtures.py +++ b/src/tespy/tools/fluid_properties/mixtures.py @@ -13,6 +13,8 @@ import math +from CoolProp.CoolProp import HAPropsSI + from tespy.tools.global_vars import FLUID_ALIASES from tespy.tools.global_vars import gas_constants @@ -36,7 +38,7 @@ def h_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): def h_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) mass_fractions_gas, molar_fraction_gas, mass_liquid, _, p_sat, pp_water = cond_check(p, T, fluid_data, water_alias) @@ -88,6 +90,56 @@ def h_mix_pT_incompressible(p, T, fluid_data, **kwargs): return h +def w_mix_fluid_data(fluid_data): + + water_alias = _get_fluid_alias("H2O", fluid_data) + water_alias = next(iter(water_alias)) + + air_alias = _get_fluid_alias("air", fluid_data) + air_alias = next(iter(air_alias)) + + return ( + fluid_data[water_alias]["mass_fraction"] + / fluid_data[air_alias]["mass_fraction"] + ) + +def w_mix_pTrh_humidair(p, T, rh): + return HAPropsSI("W", "P", p, "T", T, "RH", rh) # kg water/kg dry air + +def w_mix_pT_humidair(p, T, fluid_data, **kwargs): + w_def = w_mix_fluid_data(fluid_data) + w_max = w_mix_pTrh_humidair(p, T, 1.0) + if w_def > w_max: + _msg = f"Humidity ratio {w_def:.4f} exceeds maximum value of {w_max:.4f} for given p and T. Check fluid composition." + return w_max + return w_def + +def h_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_pT_humidair(p, T, fluid_data, **kwargs) + return HAPropsSI("H", "P", p, "T", T, "W", w) + +def w_mix_phrh_humidair(p, h, rh): + return HAPropsSI("W", "P", p, "H", h, "RH", rh) # kg water/kg dry air + +def w_mix_ph_humidair(p, h, fluid_data, **kwargs): + w_def = w_mix_fluid_data(fluid_data) + w_max = w_mix_phrh_humidair(p, h, 1.0) + if w_def > w_max: + _msg = f"Humidity ratio {w_def:.4f} exceeds maximum value of {w_max:.4f} for given p and T. Check fluid composition." + return w_max + return w_def + +def w_mix_psrh_humidair(p, s, rh): + return HAPropsSI("W", "P", p, "S", s, "RH", rh) # kg water/kg dry air + +def w_mix_ps_humidair(p, s, fluid_data, **kwargs): + w_def = w_mix_fluid_data(fluid_data) + w_max = w_mix_psrh_humidair(p, s, 1.0) + if w_def > w_max: + _msg = f"Humidity ratio {w_def:.4f} exceeds maximum value of {w_max:.4f} for given p and T. Check fluid composition." + return w_max + return w_def + def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): molar_fractions = get_molar_fractions(fluid_data) @@ -103,7 +155,7 @@ def s_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): def s_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) mass_fractions_gas, molar_fraction_gas, mass_liquid, _, p_sat, pp_water = cond_check(p, T, fluid_data, water_alias) @@ -136,6 +188,11 @@ def s_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return s +def s_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_fluid_data(fluid_data) + return HAPropsSI("S", "P", p, "T", T, "W", w) + + def v_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): molar_fractions = get_molar_fractions(fluid_data) @@ -151,7 +208,7 @@ def v_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): def v_mix_pT_ideal_cond(p=None, T=None, fluid_data=None, **kwargs): - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) _, molar_fraction_gas, mass_liquid, _, p_sat, pp_water = cond_check(p, T, fluid_data, water_alias) @@ -183,6 +240,11 @@ def v_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return v +def v_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_fluid_data(fluid_data) + return HAPropsSI("V", "P", p, "T", T, "W", w) + + def viscosity_mix_pT_ideal(p=None, T=None, fluid_data=None, **kwargs): r""" Calculate dynamic viscosity from pressure and temperature. @@ -239,10 +301,15 @@ def viscosity_mix_pT_incompressible(p=None, T=None, fluid_data=None, **kwargs): return viscosity +def viscosity_mix_pT_humidair(p, T, fluid_data, **kwargs): + w = w_mix_fluid_data(fluid_data) + return HAPropsSI("Visc", "P", p, "T", T, "W", w) + + def exergy_chemical_ideal_cond(pamb, Tamb, fluid_data, Chem_Ex): molar_fractions = get_molar_fractions(fluid_data) - water_alias = _water_in_mixture(fluid_data) + water_alias = _get_fluid_alias("H2O", fluid_data) if water_alias: water_alias = next(iter(water_alias)) _, molar_fractions_gas, _, molar_liquid, _, _ = cond_check( @@ -277,9 +344,9 @@ def exergy_chemical_ideal_cond(pamb, Tamb, fluid_data, Chem_Ex): return ex_chemical * 1e3 # Data from Chem_Ex are in kJ / mol -def _water_in_mixture(fluid_data): +def _get_fluid_alias(fluid, fluid_data): return ( - FLUID_ALIASES.get_fluid("H2O") + FLUID_ALIASES.get_fluid(fluid) & set([ f for f in fluid_data if _is_larger_than_precision(fluid_data[f]["mass_fraction"]) @@ -351,14 +418,16 @@ def cond_check(p, T, fluid_data, water_alias): T_MIX_PH_REVERSE = { "ideal": h_mix_pT_ideal, "ideal-cond": h_mix_pT_ideal_cond, - "incompressible": h_mix_pT_incompressible + "incompressible": h_mix_pT_incompressible, + "humidair": h_mix_pT_humidair } T_MIX_PS_REVERSE = { "ideal": s_mix_pT_ideal, "ideal-cond": s_mix_pT_ideal_cond, - "incompressible": s_mix_pT_incompressible + "incompressible": s_mix_pT_incompressible, + "humidair": s_mix_pT_humidair } @@ -366,28 +435,32 @@ def cond_check(p, T, fluid_data, water_alias): "ideal": h_mix_pT_ideal, "ideal-cond": h_mix_pT_ideal_cond, "incompressible": h_mix_pT_incompressible, - "forced-gas": h_mix_pT_forced_gas + "forced-gas": h_mix_pT_forced_gas, + "humidair": h_mix_pT_humidair } S_MIX_PT_DIRECT = { "ideal": s_mix_pT_ideal, "ideal-cond": s_mix_pT_ideal_cond, - "incompressible": s_mix_pT_incompressible + "incompressible": s_mix_pT_incompressible, + "humidair": s_mix_pT_humidair } V_MIX_PT_DIRECT = { "ideal": v_mix_pT_ideal, "ideal-cond": v_mix_pT_ideal_cond, - "incompressible": v_mix_pT_incompressible + "incompressible": v_mix_pT_incompressible, + "humidair": v_mix_pT_humidair } VISCOSITY_MIX_PT_DIRECT = { "ideal": viscosity_mix_pT_ideal, "ideal-cond": viscosity_mix_pT_ideal, - "incompressible": viscosity_mix_pT_incompressible + "incompressible": viscosity_mix_pT_incompressible, + "humidair": viscosity_mix_pT_humidair } EXERGY_CHEMICAL = { 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/src/tespy/tools/helpers.py b/src/tespy/tools/helpers.py index d1e78d381..1a114d8e5 100644 --- a/src/tespy/tools/helpers.py +++ b/src/tespy/tools/helpers.py @@ -16,6 +16,7 @@ from copy import deepcopy import pandas as pd +import numpy as np from tespy import __datapath__ from tespy.tools import logger @@ -475,9 +476,10 @@ def _solve_jacobian(obj, data, increment_filter, eq_num): for dependent, dx in data._vector_dependents[0].items(): f = data.func - obj._partial_derivative_fluid( - dependent, eq_num, f, dx, increment_filter, **data.func_params - ) + for dx in dx: + obj._partial_derivative_fluid( + dependent, eq_num, f, dx, increment_filter, **data.func_params + ) def _is_variable(var, increment_filter=None): @@ -699,6 +701,23 @@ def central_difference(function=None, parameter=None, delta=None, **kwargs): return (function(**upper) - function(**lower)) / (2 * delta) +def seeded_random_generator(seed_value: str | bytes | int) -> np.random.Generator: + """Generate a reproducible random number generator based on a seed value.""" + if isinstance(seed_value, str): + seed_value = seed_value.encode("utf-8") + if isinstance(seed_value, bytes): + seed_value = int.from_bytes(seed_value, "little", signed=False) % (2**32) + if not isinstance(seed_value, int): + raise ValueError("Seed value must be of type str, bytes or int.") + return np.random.default_rng(seed=seed_value) + + +def seeded_random(seed_value: str | bytes | int) -> float: + """Generate a reproducible random number between 0 and 1 based on a seed value.""" + rng = seeded_random_generator(seed_value) + return rng.random() + + def get_basic_path(): """ Return the basic tespy path and creates it if necessary. diff --git a/src/tespy/tools/optimization.py b/src/tespy/tools/optimization.py index 8150729da..5d7cbd340 100644 --- a/src/tespy/tools/optimization.py +++ b/src/tespy/tools/optimization.py @@ -293,8 +293,8 @@ def fitness(self, x): (-1) ** (sense + 1) * f for f, sense in zip(fitness, self.minimize) ] - cu = self.collect_constraints("upper") - cl = self.collect_constraints("lower") + cu = self._evaluate_constraints("upper") + cl = self._evaluate_constraints("lower") return fitness + cu + cl diff --git a/tests/test_connections.py b/tests/test_connections.py index 9fe6c5458..1a3702047 100644 --- a/tests/test_connections.py +++ b/tests/test_connections.py @@ -29,10 +29,12 @@ from tespy.connections import Ref from tespy.connections.connection import ConnectionBase from tespy.connections.connection import connection_registry +from tespy.connections.humidairconnection import HAConnection from tespy.networks import Network 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 @@ -515,7 +517,7 @@ def test_all_classes_in_registry(obj): def make_connection(cls): - if cls == Connection: + if cls == Connection or cls == HAConnection: return cls(Source(""), "out1", Sink(""), "in1") elif cls == PowerConnection: return cls(PowerSource(""), "power", PowerSink(""), "power") @@ -527,6 +529,7 @@ def make_connection(cls): QUANTITY_EXEMPTIONS = {} + def properties_of(instance): return [ prop @@ -534,6 +537,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 +547,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 +559,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/humidair.ipynb b/tests/test_tools/test_fluid_properties/humidair.ipynb new file mode 100644 index 000000000..d93eef43a --- /dev/null +++ b/tests/test_tools/test_fluid_properties/humidair.ipynb @@ -0,0 +1,187 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "5796d5ee", + "metadata": {}, + "outputs": [], + "source": [ + "from tespy.connections import HAConnection, Connection, Ref\n", + "from tespy.components import Source, Sink, MovingBoundaryHeatExchanger\n", + "from tespy.networks import Network\n", + "\n", + "from CoolProp.CoolProp import HAPropsSI" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "947c79c2", + "metadata": {}, + "outputs": [], + "source": [ + "def get_water_air_mixture_fractions_pTR(p, T, R):\n", + " w = HAPropsSI(\"W\", \"P\", p, \"T\", T, \"R\", R)\n", + " return get_water_air_mixture_fractions_w(w)\n", + "\n", + "\n", + "def get_water_air_mixture_fractions_w(w):\n", + " air = 1 / (1 + w)\n", + " return {\"Air\": air, \"water\": 1 - air}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7351d18", + "metadata": {}, + "outputs": [], + "source": [ + "from tespy.tools.fluid_properties.functions import _get_humid_air_humidity_ratio\n", + "\n", + "\n", + "def _calculate_dQ_freezing_water(T_evap, T_air, air_fluid_data, air_pressure, air_mass_flow):\n", + " dh_water = 4000 * (T_air - 273.15)\n", + " dh_freeze = 333000\n", + " dh_ice = 2000 * (273.15 - T_evap)\n", + " w_actual = _get_humid_air_humidity_ratio(air_fluid_data)\n", + " w_saturated = HAPropsSI(\"W\", \"P\", air_pressure, \"T\", T_air, \"R\", 1)\n", + "\n", + " m_water_freeze = (w_actual - w_saturated) * air_mass_flow\n", + " return (dh_water + dh_freeze + dh_ice) * m_water_freeze\n", + "\n", + "\n", + "class MovingBoundaryHeatExchangerWithFreeze(MovingBoundaryHeatExchanger):\n", + "\n", + " def energy_balance_func(self):\n", + " residual = super().energy_balance_func()\n", + " T_evap = self.outl[1].calc_T()\n", + " if T_evap < 273.15:\n", + " T_air_out = self.outl[0].calc_T()\n", + " if T_air_out > 273.15:\n", + " dQ_freeze = _calculate_dQ_freezing_water(\n", + " T_evap, T_air_out, self.outl[0].fluid_data, self.outl[0].p.val_SI, self.inl[0].m.val_SI\n", + " )\n", + " residual -= dQ_freeze\n", + "\n", + " return residual\n", + "\n", + " def energy_balance_dependents(self):\n", + " return [\n", + " self.inl[0].m, self.inl[0].h,\n", + " self.outl[0].p, self.outl[0].h,\n", + " self.inl[1].m, self.inl[1].h,\n", + " self.outl[1].p, self.outl[1].p\n", + " ]\n", + "\n", + " def calc_parameters(self):\n", + " super().calc_parameters()\n", + " self.Q.val_SI = -self.inl[1].m.val_SI * (self.outl[1].h.val_SI - self.inl[1].h.val_SI)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10a7b936", + "metadata": {}, + "outputs": [], + "source": [ + "nw = Network()\n", + "nw.units.set_defaults(\n", + " temperature=\"°C\",\n", + " pressure=\"bar\",\n", + " heat=\"kW\"\n", + ")\n", + "\n", + "\n", + "so = Source(\"source\")\n", + "hex = MovingBoundaryHeatExchanger(\"heat exchanger\")\n", + "si = Sink(\"sink\")\n", + "\n", + "so2 = Source(\"refrigerant source\")\n", + "si2 = Sink(\"refrigerant sink\")\n", + "\n", + "a1 = Connection(so2, \"out1\", hex, \"in2\", label=\"a1\")\n", + "a2 = Connection(hex, \"out2\", si2, \"in1\", label=\"a2\")\n", + "\n", + "c1 = HAConnection(so, \"out1\", hex, \"in1\", label=\"c1\")\n", + "c2 = HAConnection(hex, \"out1\", si, \"in1\", label=\"c2\")\n", + "\n", + "nw.add_conns(a1, a2, c1, c2)\n", + "\n", + "a1.set_attr(fluid={\"NH3\": 1}, m=2, x=0.3)\n", + "a2.set_attr(td_dew=5, T_dew=-20)\n", + "\n", + "c1.set_attr(\n", + " p=1,\n", + " T=15,\n", + " w0=0.004, # -> implicitly sets starting value for composition\n", + " #w=0.004, -> use this to directly impose fluid composition, then do not specify r and fluid balance for this example\n", + " r=1,\n", + " fluid_balance=True,\n", + ")\n", + "c2.set_attr(T=5)\n", + "\n", + "hex.set_attr(dp1=0, dp2=0)\n", + "\n", + "nw.solve(\"design\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3fc42bdf", + "metadata": {}, + "outputs": [], + "source": [ + "nw.print_results()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd871dc0", + "metadata": {}, + "outputs": [], + "source": [ + "heat, T_hot, T_cold, _, _ = hex.calc_sections()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f9ed087", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "\n", + "\n", + "plt.plot(heat, T_hot)\n", + "plt.plot(heat, T_cold)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "tespy", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} 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) diff --git a/tutorial/advanced/stepwise.py b/tutorial/advanced/stepwise.py index d720ffb5e..96096de4c 100644 --- a/tutorial/advanced/stepwise.py +++ b/tutorial/advanced/stepwise.py @@ -7,7 +7,7 @@ temperature="degC", pressure="bar", enthalpy="kJ/kg", power="kW", heat="kW" ) # %%[sec_2] -from tespy.components import Condenser +from tespy.components import MovingBoundaryHeatExchanger from tespy.components import CycleCloser from tespy.components import SimpleHeatExchanger from tespy.components import Pump @@ -20,7 +20,7 @@ va = Sink("valve") # consumer system -cd = Condenser("condenser") +cd = MovingBoundaryHeatExchanger("condenser") rp = Pump("recirculation pump") cons = SimpleHeatExchanger("consumer") # %%[sec_3] @@ -40,9 +40,8 @@ rp.set_attr(eta_s=0.75) cons.set_attr(pr=0.99) # %%[sec_5] -from CoolProp.CoolProp import PropsSI as PSI -p_cond = PSI("P", "Q", 1, "T", 273.15 + 95, working_fluid) / 1e5 -c0.set_attr(T=170, p=p_cond, fluid={working_fluid: 1}) +c0.set_attr(T=170, T_dew=85, fluid={working_fluid: 1}) +c1.set_attr(x=0) c20.set_attr(T=60, p=2, fluid={"water": 1}) c22.set_attr(T=90) @@ -88,11 +87,12 @@ ev.set_attr(pr1=0.99) su.set_attr(pr1=0.99, pr2=0.99) # %%[sec_10] +# condenser outlet state +c1.set_attr(x=0) # evaporator system cold side c4.set_attr(x=0.9, T=5) -h_sat = PSI("H", "Q", 1, "T", 273.15 + 15, working_fluid) / 1e3 -c6.set_attr(h=h_sat) +c6.set_attr(td_dew=5) # evaporator system hot side c17.set_attr(T=15, fluid={"water": 1}) @@ -139,10 +139,10 @@ ic.set_attr(pr1=0.99, pr2=0.98) hsp.set_attr(eta_s=0.75) # %%[sec_15] -c0.set_attr(p=p_cond, fluid={working_fluid: 1}) +c0.set_attr(T_dew=85, fluid={working_fluid: 1}) -c6.set_attr(h=c5.h.val + 10) -c8.set_attr(h=c5.h.val + 10) +c6.set_attr(td_dew=5) +c8.set_attr(td_dew=5) c7.set_attr(h=c5.h.val * 1.2) c9.set_attr(h=c5.h.val * 1.2) @@ -152,13 +152,13 @@ # %% [sec_16] nw.solve("design") -c0.set_attr(p=None) -cd.set_attr(ttd_u=5) +c0.set_attr(T_dew=None) +cd.set_attr(td_pinch=10) c4.set_attr(T=None) ev.set_attr(ttd_l=5) -c6.set_attr(h=None) +c6.set_attr(td_dew=None) su.set_attr(ttd_u=5) c7.set_attr(h=None) @@ -179,7 +179,9 @@ cons.set_attr(design=["pr"], offdesign=["zeta"]) cd.set_attr( - design=["pr2", "ttd_u"], offdesign=["zeta2", "kA_char"] + design=["pr2", "td_pinch"], offdesign=["zeta2", "UA_cecchinato"], + alpha_ratio=1, area_ratio=1, re_exp_r=0.8, re_exp_sf=0.50, + refrigerant_index=0 ) from tespy.tools.characteristics import CharLine