Steam Turbine with wet expansion

docs/whats_new/v0-7-9.rst

@@ -7,6 +7,11 @@ New Features
   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 <https://github.com/oemof/tespy/pull/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 <https://github.com/oemof/tespy/pull/602>`__).
 
 Contributors
 ############

pyproject.toml

@@ -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"}

src/tespy/components/__init__.py

@@ -26,4 +26,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

src/tespy/components/turbomachinery/steam_turbine.py

@@ -0,0 +1,294 @@
+# -*- 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.document_models import generate_latex_eq
+from tespy.tools.fluid_properties import isentropic
+from tespy.tools.fluid_properties import h_mix_pQ
+
+
+@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.
+
+    >>> 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,
+            latex=self.eta_s_wet_func_doc)
+
+        return params
+
+    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
+            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)
+
+    def eta_s_wet_func_doc(self, label):
+        r"""
+        Equation for given dry isentropic efficiency of a turbine.
+
+        Parameters
+        ----------
+        label : str
+            Label for equation.
+
+        Returns
+        -------
+        latex : str
+            LaTeX code of equations applied.
+        """
+        latex = (
+            r'\begin{split}' + '\n'
+            r'0 &=-\left(h_\mathrm{out}-h_\mathrm{in}\right)+\left('
+            r'h_\mathrm{out,s}-h_\mathrm{in}\right)\cdot\eta_\mathrm{s}\\'
+            + '\n'
+            r'\eta_\mathrm{s} &=\eta_\mathrm{s}^\mathrm{dry}\cdot \left( 1 - '
+            r'\alpha \cdot y_\mathrm{m} \right)\\' + '\n'
+            r'y_\mathrm{m} &=\frac{\left( 1-x_\mathrm{in} \right)+\left( '
+            r'1-x_\mathrm{out} \right)}{2}\\' + '\n'
+            r'\end{split}'
+        )
+        return generate_latex_eq(self, latex, label)

src/tespy/components/turbomachinery/turbine.py

@@ -12,6 +12,7 @@
 
 import numpy as np
 
+
 from tespy.components.component import component_registry
 from tespy.components.turbomachinery.base import Turbomachine
 from tespy.tools import logger
@@ -27,6 +28,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 +251,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 +526,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"""