diff --git a/docs/references.bib b/docs/references.bib index 92af8f165..3e2cada42 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -27,6 +27,15 @@ @misc{Kearney1988 year = {1988}, } +@book{Tanuma2017, + doi = {10.1016/C2014-0-03636-2}, + editor = {T. Tanuma}, + isbn = {9780081003145}, + publisher = {Elsevier}, + title = {Advances in Steam Turbines for Modern Power Plants}, + year = {2017}, +} + @TechReport{Lippke1995, author = {Lippke, Frank}, institution = {Sandia National Lab.}, @@ -411,3 +420,14 @@ @inproceedings{Fritz2024 author = {Fritz, Malte and Freißmann, Jonas and Tuschy, Ilja}, title = {Open-Source Web Dashboard zur Simulation, Analyse und Bewertung von Wärmepumpen} } + +@article{bell2015, + title = {A generalized moving-boundary algorithm to predict the heat transfer rate of counterflow heat exchangers for any phase configuration}, + journal = {Applied Thermal Engineering}, + volume = {79}, + pages = {192-201}, + year = {2015}, + issn = {1359-4311}, + doi = {https://doi.org/10.1016/j.applthermaleng.2014.12.028}, + author = {Ian H. Bell and Sylvain Quoilin and Emeline Georges and James E. Braun and Eckhard A. Groll and W. Travis Horton and Vincent Lemort}, +} diff --git a/docs/whats_new.rst b/docs/whats_new.rst index 047fddcaa..74b68deeb 100644 --- a/docs/whats_new.rst +++ b/docs/whats_new.rst @@ -3,6 +3,7 @@ What's New Discover notable new features and improvements in each release +.. include:: whats_new/v0-7-9.rst .. include:: whats_new/v0-7-8-002.rst .. include:: whats_new/v0-7-8-001.rst .. include:: whats_new/v0-7-8.rst diff --git a/docs/whats_new/v0-7-9.rst b/docs/whats_new/v0-7-9.rst new file mode 100644 index 000000000..301ac5b93 --- /dev/null +++ b/docs/whats_new/v0-7-9.rst @@ -0,0 +1,38 @@ +v0.7.9 - Under development +++++++++++++++++++++++++++ + +New Features +############ +- Implement a new property for connections to report the phase of the fluid, + i.e. :code:`"l"` for liquid, :code:`"tp"` for two-phase and :code:`"g"` for + gaseous. The phase is only reported in subcritical pressure + (`PR #592 `__). +- Implement the Baumann correlation for wet expansion in steamturbines. The + feature is available in a new component class, + :py:class:`tespy.components.turbomachinery.steam_turbine.SteamTurbine`. To + use the feature check the documentation of the component class + (`PR #602 `__). +- Implement a new component class + :py:class:`tespy.components.heat_exchangers.movingboundary.MovingBoundaryHeatExchanger`, + which allows to make specification of internal pinch temperature difference + in case of heat exchanger internal phase changes on the hot or the cold side + of the component. It splits the heat exchange automatically in sections, + where in each section the phase of the hot and the cold fluid does not change. + These sections are utilized to find the minimum pinch internally, and allow + to assing heat transfer coefficients `UA` for all individual sections as well + as the total sum + (`PR #515 `__). + +Bug Fixes +######### +- Run the postprocessing only for a converged solution. Otherwise specified + values on buses, components or connections may change to nonsense, because + the calculation is based on not-converged results of the variables + (`PR #609 `__). +- Fix the fuel cell example and make sure it is tested properly + (`PR #618 `__). + +Contributors +############ +- `@tlmerbecks `__ +- Francesco Witte (`@fwitte `__) diff --git a/pyproject.toml b/pyproject.toml index ceda28881..4df60e5e6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ exclude = ["docs/_build"] [project] name = "tespy" -version = "0.7.8.post2" +version = "0.7.9" description = "Thermal Engineering Systems in Python (TESPy)" readme = "README.rst" authors = [ @@ -51,6 +51,7 @@ dependencies = [ "matplotlib>=3.2.1", "numpy>=1.13.3", "pandas>=1.3.0", + "scipy", "tabulate>=0.8.2", ] license = {text = "MIT"} diff --git a/src/tespy/__init__.py b/src/tespy/__init__.py index 85fb30466..31307c77b 100644 --- a/src/tespy/__init__.py +++ b/src/tespy/__init__.py @@ -3,7 +3,7 @@ import os __datapath__ = os.path.join(importlib.resources.files("tespy"), "data") -__version__ = '0.7.8.post2 - Newton\'s Nature' +__version__ = '0.7.9 - Newton\'s Nature' # tespy data and connections import from . import connections # noqa: F401 diff --git a/src/tespy/components/__init__.py b/src/tespy/components/__init__.py index d459df9af..f52bd89fd 100644 --- a/src/tespy/components/__init__.py +++ b/src/tespy/components/__init__.py @@ -10,6 +10,7 @@ from .heat_exchangers.base import HeatExchanger # noqa: F401 from .heat_exchangers.condenser import Condenser # noqa: F401 from .heat_exchangers.desuperheater import Desuperheater # noqa: F401 +from .heat_exchangers.movingboundary import MovingBoundaryHeatExchanger # noqa: F401 from .heat_exchangers.parabolic_trough import ParabolicTrough # noqa: F401 from .heat_exchangers.simple import HeatExchangerSimple # noqa: F401 from .heat_exchangers.simple import SimpleHeatExchanger # noqa: F401 @@ -26,4 +27,5 @@ from .subsystem import Subsystem # noqa: F401 from .turbomachinery.compressor import Compressor # noqa: F401 from .turbomachinery.pump import Pump # noqa: F401 +from .turbomachinery.steam_turbine import SteamTurbine # noqa: F401 from .turbomachinery.turbine import Turbine # noqa: F401 diff --git a/src/tespy/components/heat_exchangers/movingboundary.py b/src/tespy/components/heat_exchangers/movingboundary.py new file mode 100644 index 000000000..a95d0db24 --- /dev/null +++ b/src/tespy/components/heat_exchangers/movingboundary.py @@ -0,0 +1,583 @@ +# -*- coding: utf-8 + +"""Module of class MovingBoundaryCondenser. + + +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/heat_exchangers/movingboundary.py + +SPDX-License-Identifier: MIT +""" +import math + +import numpy as np + +from tespy.components.heat_exchangers.base import HeatExchanger +from tespy.tools.data_containers import ComponentProperties as dc_cp +from tespy.tools.fluid_properties import T_mix_ph +from tespy.tools.fluid_properties import h_mix_pQ +from tespy.tools.fluid_properties import single_fluid +from tespy.tools.global_vars import ERR +from tespy.tools.logger import logger + + +class MovingBoundaryHeatExchanger(HeatExchanger): + r""" + Class for counter current heat exchanger with UA sections. + + The heat exchanger is internally discretized into multiple sections, which + are defined by phase changes. The component assumes, that no pressure + losses occur. In principle the implementations follows :cite:`bell2015`. + + **Mandatory Equations** + + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.energy_balance_func` + + **Optional Equations** + + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.energy_balance_hot_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.kA_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.kA_char_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.ttd_u_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.ttd_l_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.ttd_min_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.eff_cold_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.eff_hot_func` + - :py:meth:`tespy.components.heat_exchangers.base.HeatExchanger.eff_max_func` + - hot side :py:meth:`tespy.components.component.Component.pr_func` + - cold side :py:meth:`tespy.components.component.Component.pr_func` + - hot side :py:meth:`tespy.components.component.Component.zeta_func` + - cold side :py:meth:`tespy.components.component.Component.zeta_func` + - hot side :py:meth:`tespy.components.component.Component.dp_func` + - cold side :py:meth:`tespy.components.component.Component.dp_func` + - hot side :py:meth:`tespy.components.component.Component.UA_func` + - cold side :py:meth:`tespy.components.component.Component.td_pinch_func` + + Inlets/Outlets + + - in1, in2 (index 1: hot side, index 2: cold side) + - out1, out2 (index 1: hot side, index 2: cold side) + + Image + + .. image:: /api/_images/HeatExchanger.svg + :alt: flowsheet of the heat exchanger + :align: center + :class: only-light + + .. image:: /api/_images/HeatExchanger_darkmode.svg + :alt: flowsheet of the heat exchanger + :align: center + :class: only-dark + + 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. + + Q : float, dict + Heat transfer, :math:`Q/\text{W}`. + + pr1 : float, dict, :code:`"var"` + Outlet to inlet pressure ratio at hot side, :math:`pr/1`. + + pr2 : float, dict, :code:`"var"` + Outlet to inlet pressure ratio at cold side, :math:`pr/1`. + + dp1 : float, dict, :code:`"var"` + Inlet to outlet pressure delta at hot side, unit is the network's + pressure unit!. + + dp2 : float, dict, :code:`"var"` + Inlet to outlet pressure delta at cold side, unit is the network's + pressure unit!. + + zeta1 : float, dict, :code:`"var"` + Geometry independent friction coefficient at hot side, + :math:`\frac{\zeta}{D^4}/\frac{1}{\text{m}^4}`. + + zeta2 : float, dict, :code:`"var"` + Geometry independent friction coefficient at cold side, + :math:`\frac{\zeta}{D^4}/\frac{1}{\text{m}^4}`. + + ttd_l : float, dict + Lower terminal temperature difference :math:`ttd_\mathrm{l}/\text{K}`. + + ttd_u : float, dict + Upper terminal temperature difference :math:`ttd_\mathrm{u}/\text{K}`. + + ttd_min : float, dict + Minumum terminal temperature difference :math:`ttd_\mathrm{min}/\text{K}`. + + eff_cold : float, dict + Cold side heat exchanger effectiveness :math:`eff_\text{cold}/\text{1}`. + + eff_hot : float, dict + Hot side heat exchanger effectiveness :math:`eff_\text{hot}/\text{1}`. + + eff_max : float, dict + Max value of hot and cold side heat exchanger effectiveness values + :math:`eff_\text{max}/\text{1}`. + + kA : float, dict + Area independent heat transfer coefficient, + :math:`kA/\frac{\text{W}}{\text{K}}`. + + kA_char : dict + Area independent heat transfer coefficient characteristic. + + kA_char1 : tespy.tools.characteristics.CharLine, dict + Characteristic line for hot side heat transfer coefficient. + + kA_char2 : tespy.tools.characteristics.CharLine, dict + Characteristic line for cold side heat transfer coefficient. + + UA : float, dict + Sum of UA in all sections of the heat exchanger. + + td_pinch : float, dict + Value of the lowest delta T between hot side and cold side at the + different sections. + + Note + ---- + The equations only apply to counter-current heat exchangers. + + Example + ------- + Water vapor should be cooled down, condensed and then further subcooled. For + this air is heated up from 15 °C to 25 °C. + + >>> from tespy.components import Source, Sink, MovingBoundaryHeatExchanger + >>> from tespy.connections import Connection + >>> from tespy.networks import Network + >>> import numpy as np + >>> nw = Network(T_unit="C", p_unit="bar") + >>> nw.set_attr(iterinfo=False) + >>> so1 = Source("vapor source") + >>> so2 = Source("air source") + >>> cd = MovingBoundaryHeatExchanger("condenser") + >>> si1 = Sink("water sink") + >>> si2 = Sink("air sink") + >>> c1 = Connection(so1, "out1", cd, "in1", label="1") + >>> c2 = Connection(cd, "out1", si1, "in1", label="2") + >>> c11 = Connection(so2, "out1", cd, "in2", label="11") + >>> c12 = Connection(cd, "out2", si2, "in1", label="12") + >>> nw.add_conns(c1, c2, c11, c12) + + To generate good guess values, first we run the simulation with fixed + pressure on the water side. The water enters at superheated vapor state + with 15 °C superheating and leaves it with 10 °C subcooling. + + >>> c1.set_attr(fluid={"Water": 1}, p=1, Td_bp=15, m=1) + >>> c2.set_attr(p=1, Td_bp=-15) + >>> c11.set_attr(fluid={"Air": 1}, p=1, T=15) + >>> c12.set_attr(p=1, T=25) + >>> nw.solve("design") + + Now we can remove the pressure specifications on the air side and impose + the minimum pinch instead, which will determine the actual water + condensation pressure. + + >>> c1.set_attr(p=None) + >>> c2.set_attr(p=None) + >>> c12.set_attr(p=None) + >>> cd.set_attr(td_pinch=5, pr1=1, pr2=1) + >>> nw.solve("design") + >>> round(c1.p.val, 3) + 0.056 + + We can also see the temperature differences in all sections of the heat + exchanger. Since the water vapor is cooled, condensed and then subcooled, + while the air does not change phase, three sections will form: + + >>> Q_sections = cd._assign_sections() + >>> T_steps_hot, T_steps_cold = cd._get_T_at_steps(Q_sections) + >>> delta_T_between_sections = [ + ... round(T_hot - T_cold, 2) + ... for T_hot, T_cold in zip(T_steps_hot, T_steps_cold) + ... ] + >>> delta_T_between_sections + [5.0, 19.75, 10.11, 25.0] + + We can see that the lowest delta T is the first one. This is the delta T + between the hot side outlet and the cold side inlet, which can also be seen + if we have a look at the network's results. + + >>> ();nw.print_results();() # doctest: +ELLIPSIS + (...) + + If we change the subcooling degree at the water outlet, the condensation + pressure and pinch will move. + + >>> c2.set_attr(Td_bp=-5) + >>> nw.solve("design") + >>> round(c1.p.val, 3) + 0.042 + >>> Q_sections = cd._assign_sections() + >>> T_steps_hot, T_steps_cold = cd._get_T_at_steps(Q_sections) + >>> delta_T_between_sections = [ + ... round(T_hot - T_cold, 2) + ... for T_hot, T_cold in zip(T_steps_hot, T_steps_cold) + ... ] + >>> delta_T_between_sections + [9.88, 14.8, 5.0, 19.88] + + Finally, in contrast to the baseclass :code:`HeatExchanger` `kA` value, the + `UA` value takes into account the heat transfer per section and calculates + the heat transfer coefficient as the sum of all sections, while the `kA` + value only takes into account the inlet and outlet temperatures and the + total heat transfer. + + >>> round(cd.kA.val) + 173307 + >>> round(cd.UA.val) + 273449 + """ + + @staticmethod + def component(): + return 'moving boundary heat exchanger' + + def get_parameters(self): + params = super().get_parameters() + params.update({ + 'U_gas_gas': dc_cp(min_val=0), + 'U_gas_twophase': dc_cp(min_val=0), + 'U_gas_liquid': dc_cp(min_val=0), + 'U_liquid_gas': dc_cp(min_val=0), + 'U_liquid_twophase': dc_cp(min_val=0), + 'U_liquid_liquid': dc_cp(min_val=0), + 'U_twophase_gas': dc_cp(min_val=0), + 'U_twophase_twophase': dc_cp(min_val=0), + 'U_twophase_liquid': dc_cp(min_val=0), + 'A': dc_cp(min_val=0), + 'UA': dc_cp( + min_val=0, num_eq=1, func=self.UA_func, deriv=self.UA_deriv + ), + 'td_pinch': dc_cp( + min_val=0, num_eq=1, func=self.td_pinch_func, + deriv=self.td_pinch_deriv, latex=None + ) + }) + return params + + @staticmethod + def _get_h_steps(c1, c2): + """Get the steps for enthalpy for a change of state from one connection + to another + + Parameters + ---------- + c1 : tespy.connections.connection.Connection + Inlet connection. + + c2 : tespy.connections.connection.Connection + Outlet connection. + + Returns + ------- + list + Steps of enthalpy of the specified connections + """ + if c1.fluid_data != c2.fluid_data: + msg = "Both connections need to utilize the same fluid data." + raise ValueError(msg) + + if c1.p.val_SI != c2.p.val_SI: + msg = ( + "This method assumes equality of pressure for the inlet and " + "the outlet connection. The pressure values provided are not " + "equal, the results may be incorrect." + ) + # change the order of connections to have c1 as the lower enthalpy + # connection (enthalpy will be rising in the list) + if c1.h.val_SI > c2.h.val_SI: + c1, c2 = c2, c1 + + h_at_steps = [c1.h.val_SI, c2.h.val_SI] + fluid = single_fluid(c1.fluid_data) + # this should be generalized to "supports two-phase" because it does + # not work with incompressibles + is_pure_fluid = fluid is not None + + if is_pure_fluid: + try: + h_sat_gas = h_mix_pQ(c1.p.val_SI, 1, c1.fluid_data) + h_sat_liquid = h_mix_pQ(c1.p.val_SI, 0, c1.fluid_data) + except (ValueError, NotImplementedError): + return h_at_steps + + if c1.h.val_SI < h_sat_liquid: + if c2.h.val_SI > h_sat_gas + ERR: + h_at_steps = [c1.h.val_SI, h_sat_liquid, h_sat_gas, c2.h.val_SI] + elif c2.h.val_SI > h_sat_liquid + ERR: + h_at_steps = [c1.h.val_SI, h_sat_liquid, c2.h.val_SI] + + elif c1.h.val_SI < h_sat_gas - ERR: + if c2.h.val_SI > h_sat_gas + ERR: + h_at_steps = [c1.h.val_SI, h_sat_gas, c2.h.val_SI] + + return h_at_steps + + @staticmethod + def _get_Q_sections(h_at_steps, mass_flow): + """Calculate the heat exchange of every section given steps of + enthalpy and mass flow. + + Parameters + ---------- + h_at_steps : list + Enthalpy values at sections (inlet, phase change points, outlet) + mass_flow : float + Mass flow value + + Returns + ------- + float + Heat exchanged between defined steps of enthalpy. + """ + return [ + (h_at_steps[i + 1] - h_at_steps[i]) * mass_flow + for i in range(len(h_at_steps) - 1) + ] + + def _assign_sections(self): + """Assign the sections of the heat exchanger + + Returns + ------- + list + List of cumulative sum of heat exchanged defining the heat exchanger + sections. + """ + h_steps_hot = self._get_h_steps(self.inl[0], self.outl[0]) + Q_sections_hot = self._get_Q_sections(h_steps_hot, self.inl[0].m.val_SI) + Q_sections_hot = [0.0] + np.cumsum(Q_sections_hot).tolist()[:-1] + + h_steps_cold = self._get_h_steps(self.inl[1], self.outl[1]) + Q_sections_cold = self._get_Q_sections(h_steps_cold, self.inl[1].m.val_SI) + Q_sections_cold = np.cumsum(Q_sections_cold).tolist() + + return sorted(Q_sections_cold + Q_sections_hot) + + def _get_T_at_steps(self, Q_sections): + """Calculate the temperature values for the provided sections. + + Parameters + ---------- + Q_sections : list + Cumulative heat exchanged from the hot side to the colde side + defining the sections of the heat exchanger. + + Returns + ------- + tuple + Lists of cold side and hot side temperature + """ + # now put the Q_sections back on the h_steps on both sides + # Since Q_sections is defined increasing we have to start back from the + # outlet of the hot side + h_steps_hot = [self.outl[0].h.val_SI + Q / self.inl[0].m.val_SI for Q in Q_sections] + h_steps_cold = [self.inl[1].h.val_SI + Q / self.inl[1].m.val_SI for Q in Q_sections] + T_steps_hot = [ + T_mix_ph(self.inl[0].p.val_SI, h, self.inl[0].fluid_data, self.inl[0].mixing_rule) + for h in h_steps_hot + ] + T_steps_cold = [ + T_mix_ph(self.inl[1].p.val_SI, h, self.inl[1].fluid_data, self.inl[1].mixing_rule) + for h in h_steps_cold + ] + return T_steps_hot, T_steps_cold + + @staticmethod + def _calc_UA_in_sections(T_steps_hot, T_steps_cold, Q_sections): + """Calculate the UA values per section of heat exchanged. + + Parameters + ---------- + T_steps_hot : list + Temperature hot side at beginning and end of sections. + + T_steps_cold : list + Temperature cold side at beginning and end of sections. + + Q_sections : list + Heat transferred in each section. + + Returns + ------- + list + Lists of UA values per section of heat exchanged. + """ + # the temperature ranges both come with increasing values + td_at_steps = [ + T_hot - T_cold + for T_hot, T_cold in zip(T_steps_hot, T_steps_cold) + ] + + td_at_steps = [abs(td) for td in td_at_steps] + td_log_in_sections = [ + (td_at_steps[i + 1] - td_at_steps[i]) + / math.log(td_at_steps[i + 1] / td_at_steps[i]) + # round is required because tiny differences may cause + # inconsistencies due to rounding errors + if round(td_at_steps[i + 1], 6) != round(td_at_steps[i], 6) + else td_at_steps[i + 1] + for i in range(len(Q_sections)) + ] + UA_in_sections = [ + abs(Q) / td_log + for Q, td_log in zip(Q_sections, td_log_in_sections) + ] + return UA_in_sections + + def calc_UA(self): + """Calculate the sum of UA for all sections in the heat exchanger + + Returns + ------- + float + Sum of UA values of all heat exchanger sections. + """ + Q_sections = self._assign_sections() + T_steps_hot, T_steps_cold = self._get_T_at_steps(Q_sections) + Q_per_section = np.diff(Q_sections) + UA_sections = self._calc_UA_in_sections(T_steps_hot, T_steps_cold, Q_per_section) + return sum(UA_sections) + + def UA_func(self, **kwargs): + r""" + Calculate heat transfer from heat transfer coefficients for + desuperheating and condensation as well as total heat exchange area. + + Returns + ------- + residual : float + Residual value of equation. + """ + return self.UA.val - self.calc_UA() + + def UA_deriv(self, increment_filter, k): + r""" + Partial derivatives of heat transfer coefficient function. + + Parameters + ---------- + increment_filter : ndarray + Matrix for filtering non-changing variables. + + k : int + Position of derivatives in Jacobian matrix (k-th equation). + """ + f = self.UA_func + for c in self.inl + self.outl: + if self.is_variable(c.m): + self.jacobian[k, c.m.J_col] = self.numeric_deriv(f, "m", c) + if self.is_variable(c.p): + self.jacobian[k, c.p.J_col] = self.numeric_deriv(f, 'p', c) + if self.is_variable(c.h): + self.jacobian[k, c.h.J_col] = self.numeric_deriv(f, 'h', c) + + def calc_td_pinch(self): + """Calculate the pinch point temperature difference + + Returns + ------- + float + Value of the pinch point temperature difference + """ + Q_sections = self._assign_sections() + T_steps_hot, T_steps_cold = self._get_T_at_steps(Q_sections) + + td_at_steps = [ + T_hot - T_cold + for T_hot, T_cold in zip(T_steps_hot, T_steps_cold) + ] + return min(td_at_steps) + + def td_pinch_func(self): + r""" + Equation for pinch point temperature difference of a condenser. + + Returns + ------- + residual : float + Residual value of equation. + + .. math:: + + 0 = td_\text{pinch} - T_\text{sat,in,1} + + T_\left( + p_\text{in,2},\left[ + h_\text{in,2} + + \frac{\dot Q_\text{cond}}{\dot m_\text{in,2}} + \right] + \right) + """ + return self.td_pinch.val - self.calc_td_pinch() + + def td_pinch_deriv(self, increment_filter, k): + """ + Calculate partial derivates of upper terminal temperature function. + + Parameters + ---------- + increment_filter : ndarray + Matrix for filtering non-changing variables. + + k : int + Position of derivatives in Jacobian matrix (k-th equation). + """ + f = self.td_pinch_func + for c in self.inl + self.outl: + if self.is_variable(c.m, increment_filter): + self.jacobian[k, c.m.J_col] = self.numeric_deriv(f, 'm', c) + if self.is_variable(c.p, increment_filter): + self.jacobian[k, c.p.J_col] = self.numeric_deriv(f, 'p', c) + if self.is_variable(c.h, increment_filter): + self.jacobian[k, c.h.J_col] = self.numeric_deriv(f, 'h', c) + + def calc_parameters(self): + super().calc_parameters() + self.UA.val = self.calc_UA() + self.td_pinch.val = self.calc_td_pinch() + if round(self.inl[0].p.val_SI) != round(self.outl[0].p.val_SI): + msg = ( + f"The {self.__class__.__name__} instance {self.label} is " + "defined for constant pressure. The identification of the " + "heat transfer sections might be wrong in case phase changes " + "are involved in the heat transfer process." + ) + logger.warning(msg) + if round(self.inl[1].p.val_SI) != round(self.outl[1].p.val_SI): + msg = ( + f"The {self.__class__.__name__} instance {self.label} is " + "defined for constant pressure. The identification of the " + "heat transfer sections might be wrong in case phase changes " + "are involved in the heat transfer process." + ) + logger.warning(msg) diff --git a/src/tespy/components/reactors/fuel_cell.py b/src/tespy/components/reactors/fuel_cell.py index a8cdf577d..403f39b8a 100644 --- a/src/tespy/components/reactors/fuel_cell.py +++ b/src/tespy/components/reactors/fuel_cell.py @@ -132,19 +132,24 @@ class FuelCell(Component): >>> water_out = Connection(fc, 'out2', water_sink, 'in1') >>> nw.add_conns(cw_in, cw_out, oxygen_in, hydrogen_in, water_out) - The fuel cell shall produce 200kW of electrical power and 200kW of heat - with an efficiency of 0.45. The thermodynamic parameters of the input - oxygen and hydrogen are given, the mass flow rates are calculated out of - the given power output. The cooling fluid is pure water. - - >>> fc.set_attr(eta=0.45, P=-200e03, Q=-200e03, pr=0.9) - >>> cw_in.set_attr(T=25, p=1, m=1, fluid={'H2O': 1}) + The fuel cell produces 200kW of electrical power with an efficiency of 0.45. + The thermodynamic parameters of the input oxygen and hydrogen are given, + the mass flow rates are calculated out of the given power output. The + temperature of the water at the outlet should be 50 °C. The cooling fluid is + pure water and is heated up from 25 °C to 40 °C. + + >>> fc.set_attr(eta=0.45, P=-200e03, pr=0.9) + >>> cw_in.set_attr(T=25, p=1, fluid={'H2O': 1}) + >>> cw_out.set_attr(T=40) >>> oxygen_in.set_attr(T=25, p=1) >>> hydrogen_in.set_attr(T=25) + >>> water_out.set_attr(T=50) >>> nw.solve('design') - >>> P_design = fc.P.val / 1e3 - >>> round(P_design, 0) - -200.0 + >>> round(cw_in.m.val, 1) + 10.2 + >>> Q = fc.Q.val / 1e3 + >>> round(Q, 0) + -642.0 >>> round(fc.eta.val, 2) 0.45 """ diff --git a/src/tespy/components/reactors/water_electrolyzer.py b/src/tespy/components/reactors/water_electrolyzer.py index 76a83f0a9..d1299a00a 100644 --- a/src/tespy/components/reactors/water_electrolyzer.py +++ b/src/tespy/components/reactors/water_electrolyzer.py @@ -176,10 +176,10 @@ class WaterElectrolyzer(Component): >>> round(el.eta.val, 1) 0.8 >>> el_cmp.set_attr(v=None) - >>> el.set_attr(P=P_design * 0.66) + >>> el.set_attr(P=P_design * 1e6 * 0.2) >>> nw.solve('offdesign', design_path='tmp') >>> round(el.eta.val, 2) - 0.88 + 0.84 >>> shutil.rmtree('./tmp', ignore_errors=True) """ diff --git a/src/tespy/components/turbomachinery/steam_turbine.py b/src/tespy/components/turbomachinery/steam_turbine.py new file mode 100644 index 000000000..486e43d91 --- /dev/null +++ b/src/tespy/components/turbomachinery/steam_turbine.py @@ -0,0 +1,287 @@ +# -*- coding: utf-8 + +"""Module of class Turbine. + + +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/turbomachinery/turbine.py + +SPDX-License-Identifier: MIT +""" + +from scipy.optimize import brentq + +from tespy.components.component import component_registry +from tespy.components.turbomachinery.turbine import Turbine +from tespy.tools.data_containers import ComponentProperties as dc_cp +from tespy.tools.data_containers import GroupedComponentProperties as dc_gcp +from tespy.tools.fluid_properties import h_mix_pQ +from tespy.tools.fluid_properties import isentropic +from tespy.tools.fluid_properties.helpers import single_fluid +from tespy.tools.helpers import fluidalias_in_list +from tespy.tools.logger import logger + + +@component_registry +class SteamTurbine(Turbine): + r""" + Class for steam turbines with wet expansion. + + **Mandatory Equations** + + - :py:meth:`tespy.components.component.Component.fluid_func` + - :py:meth:`tespy.components.component.Component.mass_flow_func` + + **Optional Equations** + + - :py:meth:`tespy.components.component.Component.pr_func` + - :py:meth:`tespy.components.turbomachinery.base.Turbomachine.energy_balance_func` + - :py:meth:`tespy.components.turbomachinery.turbine.Turbine.eta_dry_s_func` + - :py:meth:`tespy.components.turbomachinery.turbine.Turbine.eta_s_func` + - :py:meth:`tespy.components.turbomachinery.turbine.Turbine.eta_s_char_func` + - :py:meth:`tespy.components.turbomachinery.turbine.Turbine.cone_func` + + Inlets/Outlets + + - in1 + - out1 + + Image + + .. image:: /api/_images/Turbine.svg + :alt: flowsheet of the turbine + :align: center + :class: only-light + + .. image:: /api/_images/Turbine_darkmode.svg + :alt: flowsheet of the turbine + :align: center + :class: only-dark + + 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. + + P : float, dict + Power, :math:`P/\text{W}` + + eta_s : float, dict + Isentropic efficiency, :math:`\eta_s/1` + + eta_dry_s : float, dict + Dry isentropic efficiency, :math:`\eta_s/1` + + pr : float, dict, :code:`"var"` + Outlet to inlet pressure ratio, :math:`pr/1` + + eta_s_char : tespy.tools.characteristics.CharLine, dict + Characteristic curve for isentropic efficiency, provide CharLine as + function :code:`func`. + + cone : dict + Apply Stodola's cone law (works in offdesign only). + + Example + ------- + A steam turbine expands 10 kg/s of superheated steam at 550 °C and 110 bar + to 0,5 bar at the outlet. For example, it is possible to calulate the power + output and vapour content at the outlet for a given isentropic efficiency. + The :code:`SteamTurbine` class follows the implementation of the Baumann + rule :cite:`Tanuma2017` + + >>> from tespy.components import Sink, Source, SteamTurbine + >>> from tespy.connections import Connection + >>> from tespy.networks import Network + >>> nw = Network(p_unit='bar', T_unit='C', h_unit='kJ / kg', iterinfo=False) + >>> si = Sink('sink') + >>> so = Source('source') + >>> st = SteamTurbine('steam turbine') + >>> st.component() + 'steam turbine' + >>> inc = Connection(so, 'out1', st, 'in1') + >>> outg = Connection(st, 'out1', si, 'in1') + >>> nw.add_conns(inc, outg) + + In design conditions the isentropic efficiency is specified. + >>> st.set_attr(eta_s=0.9) + >>> inc.set_attr(fluid={'water': 1}, m=10, T=250, p=20) + >>> outg.set_attr(p=0.1) + >>> nw.solve('design') + >>> round(st.P.val, 0) + -7471296.0 + >>> round(outg.x.val, 3) + 0.821 + + To capture the effect of liquid drop-out on the isentropic + efficiency, the dry turbine efficiency is specified + >>> st.set_attr(eta_s=None) + >>> st.set_attr(eta_s_dry=0.9, alpha=1.0) + >>> nw.solve('design') + >>> round(st.P.val, 0) + -7009682.0 + >>> round(outg.x.val, 3) + 0.84 + """ + + @staticmethod + def component(): + return 'steam turbine' + + def get_parameters(self): + + params = super().get_parameters() + params["alpha"] = dc_cp(min_val=0.4, max_val=2.5) + params["eta_s_dry"] = dc_cp(min_val=0.0, max_val=1.0) + params["eta_s_dry_group"] = dc_gcp( + num_eq=1, elements=["alpha", "eta_s_dry"], + func=self.eta_s_wet_func, + deriv=self.eta_s_wet_deriv + ) + + return params + + def preprocess(self, num_nw_vars): + + fluid = single_fluid(self.inl[0].fluid_data) + if fluid is None: + msg = "The SteamTurbine can only be used with pure fluids." + logger.error(msg) + raise ValueError(msg) + + if not fluidalias_in_list(fluid, ["water"]): + msg = "The SteamTurbine is intended to be used with water only." + logger.warning(msg) + + return super().preprocess(num_nw_vars) + + def eta_s_wet_func(self): + r""" + Equation for given dry isentropic efficiency of a turbine under wet + expansion. + + Returns + ------- + residual : float + Residual value of equation. + + .. math:: + + 0 = -\left( h_{out} - h_{in} \right) + + \left( h_{out,s} - h_{in} \right) \cdot \eta_{s,e} + + \eta_{s,e} = \eta_{s,e}^{dry} \cdot \left( 1 - \alpha + \cdot y_m \right) + + y_m = \frac{\left( 1-x_{in}\right)+ \left( 1-x_{out} + \right)}{2} + """ + inl = self.inl[0] + outl = self.outl[0] + + state = inl.calc_phase() + if state == "tp": # two-phase or saturated vapour + + ym = 1 - (inl.calc_x() + outl.calc_x()) / 2 # average wetness + return ( + self.calc_eta_s() + - self.eta_s_dry.val * (1 - self.alpha.val * ym) + ) + + else: # superheated vapour + if outl.calc_phase() == "g": + return self.calc_eta_s() - self.eta_s_dry.val + + dp = inl.p.val_SI - outl.p.val_SI + + # compute the pressure and enthalpy at which the expansion + # enters the vapour dome + def find_sat(frac): + psat = inl.p.val_SI - frac * dp + + # calculate enthalpy under dry expansion to psat + hout_isen = isentropic( + inl.p.val_SI, + inl.h.val_SI, + psat, + inl.fluid_data, + inl.mixing_rule, + T0=inl.T.val_SI + ) + hout = ( + inl.h.val_SI - self.eta_s_dry.val + * (inl.h.val_SI - hout_isen) + ) + + # calculate enthalpy of saturated vapour at psat + hsat = h_mix_pQ(psat, 1, inl.fluid_data) + + return hout - hsat + + frac = brentq(find_sat, 1, 0) + psat = inl.p.val_SI - frac * dp + hsat = h_mix_pQ(psat, 1, inl.fluid_data) + + # calculate the isentropic efficiency for wet expansion + ym = 1 - (1.0 + outl.calc_x()) / 2 # average wetness + eta_s = self.eta_s_dry.val * (1 - self.alpha.val * ym) + + # calculate the final outlet enthalpy + hout_isen = isentropic( + psat, + hsat, + outl.p.val_SI, + inl.fluid_data, + inl.mixing_rule, + T0=inl.T.val_SI + ) + hout = hsat - eta_s * (hsat - hout_isen) + + # return residual: outlet enthalpy = calculated outlet enthalpy + return outl.h.val_SI - hout + + def eta_s_wet_deriv(self, increment_filter, k): + r""" + Partial derivatives for dry isentropic efficiency function. + + Parameters + ---------- + increment_filter : ndarray + Matrix for filtering non-changing variables. + + k : int + Position of derivatives in Jacobian matrix (k-th equation). + """ + f = self.eta_s_wet_func + i = self.inl[0] + o = self.outl[0] + if self.is_variable(i.p, increment_filter): + self.jacobian[k, i.p.J_col] = self.numeric_deriv(f, "p", i) + if self.is_variable(o.p, increment_filter): + self.jacobian[k, o.p.J_col] = self.numeric_deriv(f, "p", o) + if self.is_variable(i.h, increment_filter): + self.jacobian[k, i.h.J_col] = self.numeric_deriv(f, "h", i) + if self.is_variable(o.h, increment_filter): + self.jacobian[k, o.h.J_col] = self.numeric_deriv(f, "h", o) diff --git a/src/tespy/components/turbomachinery/turbine.py b/src/tespy/components/turbomachinery/turbine.py index ff1b8c494..941ef86f6 100644 --- a/src/tespy/components/turbomachinery/turbine.py +++ b/src/tespy/components/turbomachinery/turbine.py @@ -27,6 +27,10 @@ class Turbine(Turbomachine): r""" Class for gas or steam turbines. + The component Turbine is the parent class for the components: + + - :py:class:`tespy.components.turbomachinery.steam_turbine.SteamTurbine` + **Mandatory Equations** - :py:meth:`tespy.components.component.Component.fluid_func` @@ -246,6 +250,22 @@ def eta_s_deriv(self, increment_filter, k): if o.h.is_var and self.it == 0: self.jacobian[k, o.h.J_col] = -1 + def calc_eta_s(self): + inl = self.inl[0] + outl = self.outl[0] + return ( + (outl.h.val_SI - inl.h.val_SI) + / (isentropic( + inl.p.val_SI, + inl.h.val_SI, + outl.p.val_SI, + inl.fluid_data, + inl.mixing_rule, + T0=inl.T.val_SI + ) - inl.h.val_SI + ) + ) + def cone_func(self): r""" Equation for stodolas cone law. @@ -505,20 +525,8 @@ def calc_parameters(self): inl = self.inl[0] outl = self.outl[0] - self.eta_s.val = ( - (outl.h.val_SI - inl.h.val_SI) - / ( - isentropic( - inl.p.val_SI, - inl.h.val_SI, - outl.p.val_SI, - inl.fluid_data, - inl.mixing_rule, - T0=inl.T.val_SI - ) - - inl.h.val_SI - ) - ) + self.eta_s.val = self.calc_eta_s() + self.pr.val = outl.p.val_SI / inl.p.val_SI def exergy_balance(self, T0): r""" diff --git a/src/tespy/connections/connection.py b/src/tespy/connections/connection.py index c80b62f9e..9e4137523 100644 --- a/src/tespy/connections/connection.py +++ b/src/tespy/connections/connection.py @@ -31,6 +31,7 @@ from tespy.tools.fluid_properties import dv_mix_pdh from tespy.tools.fluid_properties import h_mix_pQ from tespy.tools.fluid_properties import h_mix_pT +from tespy.tools.fluid_properties import phase_mix_ph from tespy.tools.fluid_properties import s_mix_ph from tespy.tools.fluid_properties import v_mix_ph from tespy.tools.fluid_properties import viscosity_mix_ph @@ -266,6 +267,7 @@ def __init__(self, source, outlet_id, target, inlet_id, if hasattr(v, "func") and v.func is not None } self.state = dc_simple() + self.phase = dc_simple() self.property_data0 = [x + '0' for x in self.property_data.keys()] self.__dict__.update(self.property_data) self.mixing_rule = None @@ -704,6 +706,7 @@ def get_parameters(self): "v_ref": dc_ref( func=self.v_ref_func, deriv=self.v_ref_deriv, num_eq=1 ), + } def build_fluid_data(self): @@ -932,6 +935,12 @@ def calc_results(self): self.x.val_SI = self.calc_x() except ValueError: self.x.val_SI = np.nan + + try: + self.phase.val = self.calc_phase() + except ValueError: + self.phase.val = np.nan + try: if not self.Td_bp.is_set: self.Td_bp.val_SI = self.calc_Td_bp() @@ -1126,6 +1135,12 @@ def get_chemical_exergy(self, pamb, Tamb, Chem_Ex): self.Ex_chemical = self.m.val_SI * self.ex_chemical + def calc_phase(self): + try: + return phase_mix_ph(self.p.val_SI, self.h.val_SI, self.fluid_data) + except NotImplementedError: + return np.nan + class Ref: r""" diff --git a/src/tespy/networks/network.py b/src/tespy/networks/network.py index b32842ca3..c57adc9d7 100644 --- a/src/tespy/networks/network.py +++ b/src/tespy/networks/network.py @@ -1109,9 +1109,11 @@ def presolve_fluid_topology(self): main_conn.fluid.val.update(fixed_fractions) main_conn.fluid.is_set = {f for f in fixed_fractions} main_conn.fluid.is_var = variable - num_var = len(variable) - for f in variable: - main_conn.fluid.val[f]: (1 - mass_fraction_sum) / num_var + # this seems to be a problem in some cases, e.g. the basic + # gas turbine tutorial + # num_var = len(variable) + # for f in variable: + # main_conn.fluid.val[f] = (1 - mass_fraction_sum) / num_var [c.build_fluid_data() for c in all_connections] for fluid in main_conn.fluid.is_var: @@ -1187,7 +1189,7 @@ def init_set_properties(self): self.all_fluids = set(self.all_fluids) cols = ( [col for prop in properties for col in [prop, f"{prop}_unit"]] - + list(self.all_fluids) + + list(self.all_fluids) + ['phase'] ) self.results['Connection'] = pd.DataFrame(columns=cols, dtype='float64') # include column for fluid balance in specs dataframe @@ -1992,8 +1994,6 @@ def solve(self, mode, init_path=None, design_path=None, logger.error(msg) return - self.postprocessing() - if not self.progress: msg = ( 'The solver does not seem to make any progress, aborting ' @@ -2005,6 +2005,8 @@ def solve(self, mode, init_path=None, design_path=None, logger.warning(msg) return + self.postprocessing() + msg = 'Calculation complete.' logger.info(msg) return @@ -2476,7 +2478,10 @@ def process_connections(self): ] + [ c.fluid.val[fluid] if fluid in c.fluid.val else np.nan for fluid in self.all_fluids + ] + [ + c.phase.val ] + ) def process_components(self): @@ -2571,7 +2576,7 @@ def print_results(self, colored=True, colors=None, print_results=True): ) # connection properties - df = self.results['Connection'].loc[:, ['m', 'p', 'h', 'T']].copy() + df = self.results['Connection'].loc[:, ['m', 'p', 'h', 'T', 'x', 'phase']].copy() df = df.astype(str) for c in df.index: if not self.get_conn(c).printout: diff --git a/src/tespy/tools/document_models.py b/src/tespy/tools/document_models.py index 67fefb37d..0db50a2f9 100644 --- a/src/tespy/tools/document_models.py +++ b/src/tespy/tools/document_models.py @@ -333,6 +333,8 @@ def document_connection_params(nw, df, specs, eqs, c, rpt): unit = col + '_unit' if col == 'Td_bp': unit = 'T_unit' + elif col == 'phase': + continue col_header = ( col.replace('_', r'\_') + ' in ' + hlp.latex_unit(nw.get_attr(unit))) diff --git a/src/tespy/tools/fluid_properties/__init__.py b/src/tespy/tools/fluid_properties/__init__.py index 987173bd5..6351f42e7 100644 --- a/src/tespy/tools/fluid_properties/__init__.py +++ b/src/tespy/tools/fluid_properties/__init__.py @@ -13,6 +13,7 @@ from .functions import h_mix_pQ # noqa: F401 from .functions import h_mix_pT # noqa: F401 from .functions import isentropic # noqa: F401 +from .functions import phase_mix_ph # noqa: F401 from .functions import s_mix_ph # noqa: F401 from .functions import s_mix_pT # noqa: F401 from .functions import v_mix_ph # noqa: F401 diff --git a/src/tespy/tools/fluid_properties/functions.py b/src/tespy/tools/fluid_properties/functions.py index f97ce42d1..0ef9d582e 100644 --- a/src/tespy/tools/fluid_properties/functions.py +++ b/src/tespy/tools/fluid_properties/functions.py @@ -155,6 +155,14 @@ def Q_mix_ph(p, h, fluid_data, mixing_rule=None): msg = "Saturation function cannot be called on mixtures." raise ValueError(msg) +def phase_mix_ph(p, h, fluid_data, mixing_rule=None): + if get_number_of_fluids(fluid_data) == 1: + pure_fluid = get_pure_fluid(fluid_data) + return pure_fluid["wrapper"].phase_ph(p, h) + else: + msg = "State function cannot be called on mixtures." + raise ValueError(msg) + def p_sat_T(T, fluid_data, mixing_rule=None): if get_number_of_fluids(fluid_data) == 1: diff --git a/src/tespy/tools/fluid_properties/wrappers.py b/src/tespy/tools/fluid_properties/wrappers.py index 6c298d54f..fd1d7274e 100644 --- a/src/tespy/tools/fluid_properties/wrappers.py +++ b/src/tespy/tools/fluid_properties/wrappers.py @@ -12,6 +12,7 @@ """ import CoolProp as CP +import numpy as np from tespy.tools.global_vars import ERR @@ -93,6 +94,9 @@ def p_sat(self, T): def Q_ph(self, p, h): self._not_implemented() + def phase_ph(self, p, h): + self._not_implemented() + def d_ph(self, p, h): self._not_implemented() @@ -227,7 +231,28 @@ def p_sat(self, T): def Q_ph(self, p, h): p = self._make_p_subcritical(p) self.AS.update(CP.HmassP_INPUTS, h, p) - return self.AS.Q() + + if self.AS.phase() == CP.iphase_twophase: + return self.AS.Q() + elif self.AS.phase() == CP.iphase_liquid: + return 0 + elif self.AS.phase() == CP.iphase_gas: + return 1 + else: # all other phases - though this should be unreachable as p is sub-critical + return -1 + + def phase_ph(self, p, h): + p = self._make_p_subcritical(p) + self.AS.update(CP.HmassP_INPUTS, h, p) + + if self.AS.phase() == CP.iphase_twophase: + return "tp" + elif self.AS.phase() == CP.iphase_liquid: + return "l" + elif self.AS.phase() == CP.iphase_gas: + return "g" + else: # all other phases - though this should be unreachable as p is sub-critical + return "state not recognised" def d_ph(self, p, h): self.AS.update(CP.HmassP_INPUTS, h, p) @@ -355,6 +380,20 @@ def Q_ph(self, p, h): p = self._make_p_subcritical(p) return self.AS(h=h / 1e3, P=p / 1e6).x + def phase_ph(self, p, h): + p = self._make_p_subcritical(p) + + phase = self.AS(h=h / 1e3, P=p / 1e6).phase + + if phase in ["Liquid"]: + return "l" + elif phase in ["Vapour"]: + return "g" + elif phase in ["Two phases", "Saturated vapor", "Saturated liquid"]: + return "tp" + else: # to ensure consistent behaviour to CoolPropWrapper + return "phase not recognised" + def d_ph(self, p, h): return self.AS(h=h / 1e3, P=p / 1e6).rho diff --git a/src/tespy/tools/global_vars.py b/src/tespy/tools/global_vars.py index 8cd773f1a..055959143 100644 --- a/src/tespy/tools/global_vars.py +++ b/src/tespy/tools/global_vars.py @@ -97,7 +97,8 @@ 'units': {'J / kgK': 1, 'kJ / kgK': 1e3, 'MJ / kgK': 1e6}, 'latex_eq': r'0 = s_\mathrm{spec} - s\left(p, h \right)', 'documentation': {'float_fmt': '{:,.2f}'} - } + }, + } combustion_gases = ['methane', 'ethane', 'propane', 'butane', 'hydrogen', 'nDodecane'] diff --git a/src/tespy/tools/helpers.py b/src/tespy/tools/helpers.py index 7b9601b43..4ea3a9b03 100644 --- a/src/tespy/tools/helpers.py +++ b/src/tespy/tools/helpers.py @@ -133,7 +133,6 @@ def convert_from_SI(property, SI_value, unit): if property == 'T': converters = fluid_property_data['T']['units'][unit] return SI_value / converters[1] - converters[0] - else: return SI_value / fluid_property_data[property]['units'][unit] diff --git a/tests/test_components/test_turbomachinery.py b/tests/test_components/test_turbomachinery.py index cb727ac84..9035a9581 100644 --- a/tests/test_components/test_turbomachinery.py +++ b/tests/test_components/test_turbomachinery.py @@ -16,6 +16,7 @@ from tespy.components import Pump from tespy.components import Sink from tespy.components import Source +from tespy.components import SteamTurbine from tespy.components import Turbine from tespy.components.turbomachinery.base import Turbomachine from tespy.connections import Connection @@ -257,12 +258,22 @@ def test_Turbine(self, tmp_path): # design value of isentropic efficiency eta_s_d = ( - (self.c2.h.val_SI - self.c1.h.val_SI) / ( - isentropic(self.c1.p.val_SI, self.c1.h.val_SI, self.c2.p.val_SI, self.c1.fluid_data, self.c1.mixing_rule) - - self.c1.h.val_SI)) - msg = ('Value of isentropic efficiency must be ' + - str(round(eta_s_d, 3)) + ', is ' + str(instance.eta_s.val) + - '.') + (self.c2.h.val_SI - self.c1.h.val_SI) + / ( + isentropic( + self.c1.p.val_SI, + self.c1.h.val_SI, + self.c2.p.val_SI, + self.c1.fluid_data, + self.c1.mixing_rule + ) + - self.c1.h.val_SI + ) + ) + msg = ( + f'Value of isentropic efficiency must be {round(eta_s_d, 3)}, is ' + f'{instance.eta_s.val}.' + ) assert round(eta_s_d, 3) == round(instance.eta_s.val, 3), msg # trigger invalid value for isentropic efficiency @@ -270,11 +281,22 @@ def test_Turbine(self, tmp_path): self.nw.solve('design') self.nw._convergence_check() eta_s = ( - (self.c2.h.val_SI - self.c1.h.val_SI) / ( - isentropic(self.c1.p.val_SI, self.c1.h.val_SI, self.c2.p.val_SI, self.c1.fluid_data, self.c1.mixing_rule) - - self.c1.h.val_SI)) - msg = ('Value of isentropic efficiency must be ' + str(eta_s) + - ', is ' + str(instance.eta_s.val) + '.') + (self.c2.h.val_SI - self.c1.h.val_SI) + / ( + isentropic( + self.c1.p.val_SI, + self.c1.h.val_SI, + self.c2.p.val_SI, + self.c1.fluid_data, + self.c1.mixing_rule + ) + - self.c1.h.val_SI + ) + ) + msg = ( + f'Value of isentropic efficiency must be {round(eta_s, 3)}, is ' + f'{instance.eta_s.val}.' + ) assert round(eta_s, 3) == round(instance.eta_s.val, 3), msg # unset isentropic efficiency and inlet pressure, @@ -287,21 +309,23 @@ def test_Turbine(self, tmp_path): self.nw.solve('offdesign', design_path=tmp_path) self.nw._convergence_check() # check efficiency - msg = ('Value of isentropic efficiency (' + str(instance.eta_s.val) + - ') must be identical to design case (' + str(eta_s_d) + ').') + msg = ( + f'Value of isentropic efficiency ({instance.eta_s.val}) must be ' + f'identical to design case ({eta_s_d}).' + ) assert round(eta_s_d, 2) == round(instance.eta_s.val, 2), msg # check pressure - msg = ('Value of inlet pressure (' + str(round(self.c1.p.val_SI)) + - ') must be identical to design case (' + - str(round(self.c1.p.design)) + ').') + msg = ( + f'Value of inlet pressure ({round(self.c1.p.val_SI)}) must be ' + f'identical to design case ({round(self.c1.p.design)}).' + ) assert round(self.c1.p.design) == round(self.c1.p.val_SI), msg # lowering mass flow, inlet pressure must sink according to cone law self.c1.set_attr(m=self.c1.m.val * 0.8) self.nw.solve('offdesign', design_path=tmp_path) self.nw._convergence_check() - msg = ('Value of pressure ratio (' + str(instance.pr.val) + - ') must be at (' + str(0.128) + ').') + msg = f'Value of pressure ratio ({instance.pr.val}) must be at 0.128.' assert 0.128 == round(instance.pr.val, 3), msg # testing more parameters for eta_s_char @@ -312,28 +336,76 @@ def test_Turbine(self, tmp_path): self.nw._convergence_check() expr = self.c1.v.val_SI / self.c1.v.design eta_s = round(eta_s_d * instance.eta_s_char.char_func.evaluate(expr), 3) - msg = ('Value of isentropic efficiency (' + - str(round(instance.eta_s.val, 3)) + - ') must be (' + str(eta_s) + ').') + msg = ( + f'Value of isentropic efficiency ({round(instance.eta_s.val, 3)}) ' + f'must be {eta_s}.' + ) assert eta_s == round(instance.eta_s.val, 3), msg # test parameter specification pr instance.eta_s_char.param = 'pr' self.nw.solve('offdesign', design_path=tmp_path) self.nw._convergence_check() - expr = (self.c2.p.val_SI * self.c1.p.design / - (self.c2.p.design * self.c1.p.val_SI)) + expr = ( + self.c2.p.val_SI * self.c1.p.design + / (self.c2.p.design * self.c1.p.val_SI) + ) eta_s = round(eta_s_d * instance.eta_s_char.char_func.evaluate(expr), 3) - msg = ('Value of isentropic efficiency (' + - str(round(instance.eta_s.val, 3)) + - ') must be (' + str(eta_s) + ').') + msg = ( + f'Value of isentropic efficiency ({round(instance.eta_s.val, 3)}) ' + f'must be {eta_s}.' + ) assert eta_s == round(instance.eta_s.val, 3), msg + def test_SteamTurbine(self, tmp_path): + instance = SteamTurbine('turbine') + self.setup_network(instance) + fl = {'H2O': 1} + # start in gas, end in gas + self.c1.set_attr(fluid=fl, m=15, p=100, T=500) + self.c2.set_attr(Td_bp=1) + + eta_s_dry = 0.9 + instance.set_attr(eta_s_dry=0.9, alpha=1) + + self.nw.solve('design') + self.nw._convergence_check() + + eta_s = round(instance.eta_s.val, 4) + msg = ( + f"In gas expansion isentropic ({eta_s}) and dry " + f"efficiency ({eta_s_dry}) have to be identical." + ) + assert eta_s_dry == eta_s, msg + + # end in two phase + self.c2.set_attr(Td_bp=None, x=0.9) + self.nw.solve('design') + self.nw._convergence_check() + + eta_s = round(instance.eta_s.val, 4) + msg = ( + f"In gas expansion isentropic ({eta_s}) and dry " + f"efficiency ({eta_s_dry}) cannot be identical." + ) + assert eta_s_dry != eta_s, msg + + # start in two phase + self.c1.set_attr(T=None, x=0.99) + self.nw.solve('design') + self.nw._convergence_check() + + eta_s = round(instance.eta_s.val, 4) + msg = ( + f"In gas expansion isentropic ({eta_s}) and dry " + f"efficiency ({eta_s_dry}) cannot be identical." + ) + assert eta_s_dry != eta_s, msg + def test_Turbomachine(self): """Test component properties of turbomachines.""" instance = Turbomachine('turbomachine') - msg = ('Component name must be turbomachine, is ' + - instance.component() + '.') + msg = f'Component name must be turbomachine, is {instance.component()}.' assert 'turbomachine' == instance.component(), msg self.setup_network(instance) fl = {'N2': 0.7556, 'O2': 0.2315, 'Ar': 0.0129}