Combined cycle power plant

In this model a gas turbine with heat recovery steam generator and a steam cycle is modeled. The steam cycle has an extraction turbine, which allows to extract steam partially to provide district heating at the loss of electricity generation. The model is the combined cycle example of the exerpy library. It uses the ModelTemplate class.

On top of the model itself, the ModelTemplate is used to

• draw a Ts diagram of the steam cycle,

• draw the combined QT diagram of the heat recovery steam generator,

• run sensitivity analyses on the design parameters, and

• optimize the live steam pressure and the superheater sizing

Imports

Next to the tespy imports we need the FluidPropertyDiagram from fluprodia for the Ts diagram and the PSO (particle swarm optimization) algorithm from pymoo.

import os
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from fluprodia import FluidPropertyDiagram
from pymoo.algorithms.soo.nonconvex.pso import PSO

from tespy.components import (
    Compressor, Condenser, CycleCloser, DiabaticCombustionChamber, Drum,
    Generator, HeatExchanger, HeatSink, Merge, Motor, PowerBus, PowerSink,
    Pump, SimpleHeatExchanger, Sink, Source, Splitter, Turbine, Valve
)
from tespy.connections import Connection, HeatConnection, PowerConnection, Ref
from tespy.models.template import ModelTemplate
from tespy.networks import Network

The CombinedCyclePlant class

The model class inherits from ModelTemplate, which implements input and output parameter handling, solving, plotting diagrams and running sensitivity analyses and optimization.

We need to implement these three methods:

• _create_network to build the topology and apply the initial boundary conditions specifications.

• _parameter_lookup to map names onto parameters in the network.

• solve_model to let the optimizer and sensitivity analysis know, how to solve a model.

In _create_network the specification is applied in stages: some parameters are set, the model is solved, and then they are traded for the ones we actually want. The target specification does not converge from the default starting values, so we build the model with specifications the solver finds easy and swap them afterwards.

Next to parameters of connections and components, _parameter_lookup also provides four derived results: the fuel input from the lower heating value, the net electrical efficiency, the fuel utilisation including district heat and the stack temperature.

create_diagram is overridden because the base class draws isenthalps, which a Ts diagram of a steam cycle does not need, and the label of the zero enthalpy line cannot be placed for water.

class CombinedCyclePlant(ModelTemplate):
    """Combined cycle power plant with district heating extraction."""

    def _create_network(self):
        self.nw = nw = Network(iterinfo=False)
        nw.units.set_defaults(
            temperature="degC",
            pressure="bar",
            pressure_difference="kPa",
            enthalpy="kJ / kg",
            mass_flow="kg / s",
            heat="MW",
            power="MW",
        )

        air_in = Source("air inlet")
        fuel_in = Source("fuel inlet")
        flue_gas_out = Sink("flue gas outlet")
        compressor = Compressor("COMP")
        combustion = DiabaticCombustionChamber("CC")
        gasturbine = Turbine("GT")

        superheater = HeatExchanger("SH")
        evaporator = HeatExchanger("EVA")
        economizer = HeatExchanger("ECO")

        drum = Drum("drum")
        feed_pump = Pump("PUMP2")
        condensate_pump = Pump("PUMP1")
        drum_pump = Pump("drum pump")
        hp_steam_turbine = Turbine("ST1")
        lp_steam_turbine = Turbine("ST2")

        deaerator = Merge("DEA", num_in=3)
        dea_steam_valve = Valve("dea steam inlet valve")
        extraction = Splitter("turbine outlet extraction", num_out=3)
        heating_condenser = SimpleHeatExchanger("HC")
        main_condenser = Condenser("COND")

        water_in = Source("cooling water source")
        water_out = Sink("cooling water sink")

        cc = CycleCloser("cycle closer")

        c1 = Connection(air_in, "out1", compressor, "in1", label="1")
        c2 = Connection(compressor, "out1", combustion, "in1", label="2")
        c3 = Connection(fuel_in, "out1", combustion, "in2", label="3")
        c4 = Connection(combustion, "out1", gasturbine, "in1", label="4")
        c5 = Connection(gasturbine, "out1", superheater, "in1", label="5")
        c6 = Connection(superheater, "out1", evaporator, "in1", label="6")
        c7 = Connection(evaporator, "out1", economizer, "in1", label="7")
        c8 = Connection(economizer, "out1", flue_gas_out, "in1", label="8")

        nw.add_conns(c1, c2, c3, c4, c5, c6, c7, c8)

        c9 = Connection(superheater, "out2", cc, "in1", label="9")
        c0 = Connection(cc, "out1", hp_steam_turbine, "in1", label="0")
        c101 = Connection(hp_steam_turbine, "out1", extraction, "in1", label="101")
        c10 = Connection(extraction, "out1", dea_steam_valve, "in1", label="10")
        c10a = Connection(dea_steam_valve, "out1", deaerator, "in1", label="10a")

        c11 = Connection(extraction, "out2", heating_condenser, "in1", label="11")
        c12 = Connection(extraction, "out3", lp_steam_turbine, "in1", label="12")
        c13 = Connection(lp_steam_turbine, "out1", main_condenser, "in1", label="13")

        c14 = Connection(water_in, "out1", main_condenser, "in2", label="14")
        c15 = Connection(main_condenser, "out2", water_out, "in1", label="15")

        c16 = Connection(main_condenser, "out1", condensate_pump, "in1", label="16")

        c17 = Connection(condensate_pump, "out1", deaerator, "in3", label="17")
        c18 = Connection(heating_condenser, "out1", deaerator, "in2", label="18")
        c20 = Connection(deaerator, "out1", feed_pump, "in1", label="20")

        c21 = Connection(feed_pump, "out1", economizer, "in2", label="21")
        c22 = Connection(economizer, "out2", drum, "in1", label="22")
        c22a = Connection(drum, "out1", drum_pump, "in1", label="22a")
        c22b = Connection(drum_pump, "out1", evaporator, "in2", label="22b")
        c22c = Connection(evaporator, "out2", drum, "in2", label="22c")
        c23 = Connection(drum, "out2", superheater, "in2", label="23")

        nw.add_conns(c9, c0, c101, c10, c10a, c11, c12, c13, c14, c15, c16, c17, c18, c20, c21, c22, c22a, c22b, c22c, c23)

        gt_shaft = PowerBus("GT shaft", num_in=1, num_out=2)
        st_shaft = PowerBus("ST shaft", num_in=2, num_out=1)
        gt_generator = Generator("GEN1")
        st_generator = Generator("GEN2")
        fp_motor = Motor("MOT2")
        cp_motor = Motor("MOT1")
        dp_motor = Motor("drum pump motor")
        distribution = PowerBus("electricity distribution", num_in=2, num_out=4)
        grid = PowerSink("grid")

        e1 = PowerConnection(gasturbine, "power", gt_shaft, "power_in1", label="e01")
        e2 = PowerConnection(gt_shaft, "power_out1", compressor, "power", label="e02")
        e3 = PowerConnection(gt_shaft, "power_out2", gt_generator, "power_in", label="e03")
        e4 = PowerConnection(gt_generator, "power_out", distribution, "power_in1", label="e04")

        e5 = PowerConnection(hp_steam_turbine, "power", st_shaft, "power_in1", label="e05")
        e6 = PowerConnection(lp_steam_turbine, "power", st_shaft, "power_in2", label="e06")
        e7 = PowerConnection(st_shaft, "power_out1", st_generator, "power_in", label="e07")
        e8 = PowerConnection(st_generator, "power_out", distribution, "power_in2", label="e08")

        e9 = PowerConnection(distribution, "power_out1", fp_motor, "power_in", label="e09")
        e10 = PowerConnection(fp_motor, "power_out", feed_pump, "power", label="e10")
        e11 = PowerConnection(distribution, "power_out2", cp_motor, "power_in", label="e11")
        e12 = PowerConnection(cp_motor, "power_out", condensate_pump, "power", label="e12")
        e13 = PowerConnection(distribution, "power_out3", dp_motor, "power_in", label="e13")
        e14 = PowerConnection(dp_motor, "power_out", drum_pump, "power", label="e14")
        e15 = PowerConnection(distribution, "power_out4", grid, "power", label="e15")

        nw.add_conns(e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15)

        heating = HeatSink("heating")
        h1 = HeatConnection(heating_condenser, "heat", heating, "heat", label="h1")
        nw.add_conns(h1)

        gt_generator.set_attr(eta=0.985)
        st_generator.set_attr(eta=0.985)
        fp_motor.set_attr(eta=0.985)
        cp_motor.set_attr(eta=0.985)
        dp_motor.set_attr(eta=0.985)

        compressor.set_attr(eta_s=0.9)
        combustion.set_attr(eta=1)
        gasturbine.set_attr(eta_s=0.92)

        hp_steam_turbine.set_attr(eta_s=0.93)
        lp_steam_turbine.set_attr(eta_s=0.89)

        feed_pump.set_attr(eta_s=0.8)
        condensate_pump.set_attr(eta_s=0.8)

        evaporator.set_attr(dp1=0.7, dp2=3, ttd_l=10)
        superheater.set_attr(dp1=0.7, dp2=5)
        economizer.set_attr(dp1=0.7, dp2=2)
        drum_pump.set_attr(eta_s=0.8)
        heating_condenser.set_attr(Q=-100)
        main_condenser.set_attr(dp1=0, dp2=0)

        c1.set_attr(fluid={"AR": 0.01282, "CO2": 0.00040, "H2O": 0.00634, "N2": 0.75051, "O2": 0.22993}, p=1.013, T=15)
        c2.set_attr(p=15.51)
        c3.set_attr(fluid={"CH4": 1}, p=Ref(c2, 1, 0), T=15)
        c4.set_attr(p=15, T=1150)

        c8.set_attr(p=1.013)
        c9.set_attr(fluid={"water": 1}, p=50, T=500)
        c10.set_attr(p=15)
        c10a.set_attr(p=10)
        c13.set_attr(p=0.05)
        c14.set_attr(fluid={"water": 1}, p=1.013, T=15)
        c15.set_attr(T=28)

        c18.set_attr(td_bubble=10)
        c20.set_attr(x=0)
        c22.set_attr(td_bubble=6)
        c22c.set_attr(x=0.15)

        e15.set_attr(E=300)

        self.nw.solve("design")
        self.nw.assert_convergence()

        # a fixed live steam temperature reaches a solution from the default
        # starting values, the superheater terminal temperature difference
        # only does so from a converged state
        c9.set_attr(T=None)
        superheater.set_attr(ttd_u=25)
        self.nw.solve("design")
        self.nw.assert_convergence()
        self._stable_solution = self.nw.save(as_dict=True)

    def _parameter_lookup(self):
        return {
            "live steam pressure": ["Connections", "9", "p"],
            "compressor pressure ratio": ["Components", "COMP", "pr"],
            "air mass flow": ["Connections", "1", "m"],
            "turbine inlet temperature": ["Connections", "4", "T"],
            "condenser pressure": ["Connections", "13", "p"],
            "superheater ttd": ["Components", "SH", "ttd_u"],
            "evaporator pinch": ["Components", "EVA", "ttd_l"],
            "electricity output": ["Connections", "e15", "E"],
            "district heat": ["Components", "HC", "Q"],
            "gas turbine eta_s": ["Components", "GT", "eta_s"],
            "fuel input": {"get": self._calc_fuel_input},
            "net efficiency": {"get": self._calc_net_efficiency},
            "fuel utilisation": {"get": self._calc_fuel_utilisation},
            "stack temperature": {"get": self._calc_stack_temperature},
        }

    def _calc_fuel_input(self):
        """Thermal input from the fuel lower heating value in MW."""
        if not self._solved:
            return np.nan
        return self.nw.get_comp("CC").ti.val

    def _calc_net_efficiency(self):
        if not self._solved:
            return np.nan
        return self.nw.get_conn("e15").E.val / self.nw.get_comp("CC").ti.val

    def _calc_fuel_utilisation(self):
        """Electricity plus district heat over fuel input."""
        if not self._solved:
            return np.nan
        output = (
            self.nw.get_conn("e15").E.val
            + abs(self.nw.get_conn("h1").E.val)
        )
        return output / self.nw.get_comp("CC").ti.val

    def _calc_stack_temperature(self):
        """Flue gas temperature leaving the economizer in °C."""
        if not self._solved:
            return np.nan
        return self.nw.get_conn("8").T.val

    def create_diagram(self, fluid_name):
        """Isolines suited to a steam cycle Ts diagram.

        The base class draws isenthalps, which a Ts diagram of a steam cycle
        does not need - and the label placement of the zero enthalpy line
        cannot be resolved for water, which breaks the figure. Isobars and
        vapour quality lines are what this diagram is read for.
        """
        diagram = FluidPropertyDiagram(fluid_name)
        diagram.set_unit_system(self.nw.units)
        diagram.set_isolines_subcritical(10, 560)
        diagram.set_isolines(
            p=np.array([0.05, 0.2, 1, 5, 10, 20, 50, 100]),
            h=np.array([]),
            vol=np.array([]),
            Q=np.linspace(0, 1, 11),
        )
        diagram.calc_isolines()
        return diagram

    def solve_model(self, **kwargs):
        self.solve_model_design(**kwargs)

Design simulation

Instantiating the class builds the network and runs the staged solve, solve_model re-solves from the converged state. The electricity output is imposed at 300 MW and the district heat at 100 MW.

model = CombinedCyclePlant()
model.solve_model()

pd.Series(model.get_results([
    "fuel input", "electricity output", "district heat",
    "net efficiency", "fuel utilisation", "stack temperature",
    "live steam pressure", "superheater ttd", "evaporator pinch",
]), name="design point").round(4)

fuel input             618.6892
electricity output     300.0000
district heat         -100.0000
net efficiency           0.4849
fuel utilisation         0.6465
stack temperature      233.1641
live steam pressure     50.0000
superheater ttd         25.0000
evaporator pinch        10.0000
Name: design point, dtype: float64

Just under half of the fuel is converted to electricity. Counting the district heat as well, about two thirds of the fuel input leave the plant as useful energy. The remaining energy is lost through the chimney or in the condenser.

Ts diagram

fig, ax = model.plot_Ts_diagram_matplotlib("0")

Heat recovery steam generator

We can merge the QT diagrams of the different heat exchangers that form the HRSG. Due to the drum’s approach point temperature difference at the liquid inlet we see a jump on the water side temperature.

def hrsg_profile(model):
    """Combined flue gas and water temperature profile across the HRSG.

    Returns cumulative duty in MW together with the gas and water temperatures,
    ordered from the cold end. The saturation plateau of the evaporator is
    inserted so the pinch appears where it physically occurs.
    """
    nw = model.nw
    Q_eco, Q_eva, Q_sh = (
        abs(nw.get_comp(c).Q.val) for c in ("ECO", "EVA", "SH")
    )
    T_sat = nw.get_conn("23").T.val  # saturated steam leaving the drum

    Q = [0, Q_eco, Q_eco, Q_eco + Q_eva, Q_eco + Q_eva + Q_sh]
    gas = [nw.get_conn(c).T.val for c in ("8", "7", "7", "6", "5")]
    water = [
        nw.get_conn("21").T.val, nw.get_conn("22").T.val,
        T_sat, T_sat, nw.get_conn("9").T.val
    ]
    return np.array(Q), np.array(gas), np.array(water)


Q, gas, water = hrsg_profile(model)

fig, ax = plt.subplots(figsize=(10.5, 5.4))
ax.plot(Q, gas, "o-", color="tab:red", label="flue gas")
ax.plot(Q, water, "o-", color="tab:blue", label="water and steam")
ax.fill_between(Q, water, gas, color="tab:gray", alpha=.15)

pinch = int(np.argmin(gas - water))
ax.annotate(
    f"pinch {gas[pinch] - water[pinch]:.1f} K",
    xy=(Q[pinch], (gas[pinch] + water[pinch]) / 2),
    xytext=(12, -34), textcoords="offset points",
    arrowprops=dict(arrowstyle="->", color="gray"),
)
for x, name in zip(
        [Q[1] / 2, (Q[1] + Q[3]) / 2, (Q[3] + Q[4]) / 2],
        ["economizer", "evaporator", "superheater"]
    ):
    ax.annotate(name, xy=(x, 150), ha="center", fontsize=9, color="gray")
ax.set_xlabel("cumulative transferred heat in MW")
ax.set_ylabel("temperature in °C")
ax.legend()
ax.grid(alpha=.3)
plt.close(fig)
fig

Sensitivity analysis

At the example of the parameters below a sensitivity analysis is carried out. All simulations run in design mode, so every point is a plant designed for that specification.

parameter	design value	range
live steam pressure	50 bar	30 - 110 bar
evaporator pinch	10 K	5 - 25 K
superheater ttd	25 K	10 - 50 K

n = 6 if os.getenv("TESPY_DOCS_BUILD") else 11
design = {'live steam pressure': 50, 'superheater ttd': 25, 'evaporator pinch': 10}

pressure = model.sensitivity_analysis(
    param_dict={"live steam pressure": np.linspace(30, 110, n)},
    result_param_list=["net efficiency", "stack temperature", "fuel input"],
)
model.solve_model(**design)

pinch = model.sensitivity_analysis(
    param_dict={"evaporator pinch": np.linspace(5, 25, n)},
    result_param_list=["net efficiency", "stack temperature"],
)
model.solve_model(**design)

ttd = model.sensitivity_analysis(
    param_dict={"superheater ttd": np.linspace(10, 50, n)},
    result_param_list=["net efficiency", "stack temperature"],
)
model.solve_model(**design)

The live steam pressure has an optimum. A higher pressure improves the steam cycle, but it also raises the evaporation temperature, which leads to a higher exhaust temperature of the flue gases leaving the HRSG.

Evaporator pinch and superheater terminal temperature difference are monotonic, a lower value is (as expected) always better.

fig, ax = plt.subplots(1, 3, figsize=(13, 3.8), sharey=True)
ax[0].plot(pressure["live steam pressure"], pressure["net efficiency"], "o-", color="tab:blue")
best = pressure.loc[pressure["net efficiency"].idxmax()]
ax[0].plot(best["live steam pressure"], best["net efficiency"], "*", color="tab:red", markersize=16)
ax[0].set_xlabel("live steam pressure in bar")
ax[1].plot(pinch["evaporator pinch"], pinch["net efficiency"], "o-", color="tab:green")
ax[1].set_xlabel("evaporator pinch in K")
ax[2].plot(ttd["superheater ttd"], ttd["net efficiency"], "o-", color="tab:orange")
ax[2].set_xlabel("superheater ttd in K")
ax[0].set_ylabel("net electrical efficiency")
for a in ax:
    a.grid(alpha=.3)
plt.tight_layout()
plt.close(fig)
fig

We can check the temperature of the flue gas leaving the HRSG. The higher it is the more energy is lost to the ambient.

fig, ax = plt.subplots(figsize=(10.5, 5.2))
ax.plot(pressure["live steam pressure"], pressure["stack temperature"], "o-", color="tab:red")
ax.set_xlabel("live steam pressure in bar")
ax.set_ylabel("stack temperature in °C")
ax.grid(alpha=.3)
plt.close(fig)
fig