Supercritical CO2 recompression cycle¶
This is a model of a supercritical CO2 Brayton cycle with recompression and two recuperators providing 100 MW net power. The model is based on the model used validate the exergy analysis originally published in [9]. It is validated against the data of Penkuhn and Tsatsaronis [10]. The original model is available in the sCO2_exergy repository.
On top of the model validating the data reported, the ModelTemplate is used to implement
a log(p)-h diagram of the process,
run sensitivity analyses on some parameters, and
apply a multi-objective optimization of efficiency and recuperator UA
Imports¶
Next to the tespy imports we need the FluidPropertyDiagram from
fluprodia and the NSGA2 algorithm for
multi-objective optimization 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.moo.nsga2 import NSGA2
from tespy.components import (
Compressor, CycleCloser, Generator, HeatSource, Merge, Motor,
PowerBus, PowerSink, SectionedHeatExchanger, SimpleHeatExchanger,
Splitter, Turbine
)
from tespy.connections import Connection, HeatConnection, PowerConnection, Ref
from tespy.models.template import ModelTemplate
from tespy.networks import Network
The SCO2Model 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_networkto build the topology and apply the initial boundary conditions specifications._parameter_lookupto map names onto parameters in the network.solve_modelto let the optimizer and sensitivity analysis know, how to solve a model.
create_diagram is overridden because the base class isolines do not go beyond
200 °C while this cycle reaches 600 °C at 250 bar.
class SCO2Model(ModelTemplate):
"""Supercritical CO2 power cycle with recompression and recuperators."""
# specification of ambient state
pamb = 1.01325 # bar
Tamb = 15 # °C
def _create_network(self):
self.nw = Network(iterinfo=False)
self.nw.units.set_defaults(
temperature="°C",
pressure="bar",
pressure_difference="bar",
enthalpy="kJ / kg",
entropy="kJ / kgK",
power="MW",
heat="MW",
heat_transfer_coefficient="MW/K"
)
# components definition
closer = CycleCloser("Cycle closer")
cp1 = Compressor("Compressor 1")
cp2 = Compressor("Compressor 2")
rec1 = SectionedHeatExchanger("Recuperator 1", num_sections=10)
rec2 = SectionedHeatExchanger("Recuperator 2", num_sections=10)
cooler = SimpleHeatExchanger("Water cooler")
heater = SimpleHeatExchanger("Heater")
turb = Turbine("Turbine")
sp1 = Splitter("Splitter 1")
m1 = Merge("Merge 1")
# connections definition
# power cycle
c1 = Connection(cooler, "out1", cp1, "in1", label="01")
c2 = Connection(cp1, "out1", rec1, "in2", label="02")
c3 = Connection(rec2, "out2", heater, "in1", label="03")
c0 = Connection(heater, "out1", closer, "in1", label="00")
c4 = Connection(closer, "out1", turb, "in1", label="04")
c5 = Connection(turb, "out1", rec2, "in1", label="05")
c6 = Connection(sp1, "out1", cooler, "in1", label="06")
c10 = Connection(sp1, "out2", cp2, "in1", label="10")
c11 = Connection(cp2, "out1", m1, "in2", label="11")
c12 = Connection(rec1, "out2", m1, "in1", label="12")
c13 = Connection(m1, "out1", rec2, "in2", label="13")
c14 = Connection(rec2, "out1", rec1, "in1", label="14")
c15 = Connection(rec1, "out1", sp1, "in1", label="15")
# add connections to network
self.nw.add_conns(c0, c1, c2, c3, c4, c5, c6, c10, c11, c12, c13, c14, c15)
# power system
grid = PowerSink("grid")
distribution = PowerBus("distribution", num_in=1, num_out=3)
motor1 = Motor("motor 1")
motor2 = Motor("motor 2")
generator = Generator("generator")
e1 = PowerConnection(distribution, "power_out3", grid, "power", label="e1")
e2 = PowerConnection(distribution, "power_out1", motor1, "power_in", label="e2")
e3 = PowerConnection(motor1, "power_out", cp1, "power", label="e3")
e4 = PowerConnection(distribution, "power_out2", motor2, "power_in", label="e4")
e5 = PowerConnection(motor2, "power_out", cp2, "power", label="e5")
e6 = PowerConnection(turb, "power", generator, "power_in", label="e6")
e7 = PowerConnection(generator, "power_out", distribution, "power_in1", label="e7")
self.nw.add_conns(e1, e2, e3, e4, e5, e6, e7)
# heat input
heatsource = HeatSource("heatsource")
h1 = HeatConnection(heatsource, "heat", heater, "heat", label="h1")
self.nw.add_conns(h1)
motor1.set_attr(eta=0.97 * 0.98)
motor2.set_attr(eta=0.97 * 0.98)
generator.set_attr(eta=0.99 * 0.99)
# connection parameters
c1.set_attr(T=35, p=75, fluid={"CO2": 1})
c2.set_attr(p=258.4)
c3.set_attr(p=257)
c4.set_attr(T=600, p=250)
c5.set_attr(p=77.95)
c6.set_attr(p=75.15)
c11.set_attr(p=257.51, T=Ref(c12, 1, 0))
c14.set_attr(p=76.94)
# component parameters
turb.set_attr(eta_s=0.9)
cp1.set_attr(eta_s=0.85)
cp2.set_attr(eta_s=0.85)
rec1.set_attr(ttd_l=5)
rec2.set_attr(ttd_l=5)
# net power output
e1.set_attr(E=100)
def _parameter_lookup(self):
return {
"net power": ["Connections", "e1", "E"],
"heat input": ["Connections", "h1", "E"],
"turbine inlet temperature": ["Connections", "04", "T"],
"turbine inlet pressure": ["Connections", "04", "p"],
"compressor inlet temperature": ["Connections", "01", "T"],
"compressor inlet pressure": ["Connections", "01", "p"],
"turbine efficiency": ["Components", "Turbine", "eta_s"],
"recuperator 1 pinch": ["Components", "Recuperator 1", "td_pinch"],
"recuperator 2 pinch": ["Components", "Recuperator 2", "td_pinch"],
"recuperator 1 ttd_l": ["Components", "Recuperator 1", "ttd_l"],
"recuperator 2 ttd_l": ["Components", "Recuperator 2", "ttd_l"],
"recuperator 1 UA": ["Components", "Recuperator 1", "UA"],
"recuperator 2 UA": ["Components", "Recuperator 2", "UA"],
"recuperator UA": {"get": self._get_recuperator_total_UA},
"T recompressed in": ["Connections", "11", "T"],
"T recuperated in": ["Connections", "12", "T"],
"thermal efficiency": {"get": self._calc_thermal_efficiency},
"split ratio": {"set": self._set_split_ratio, "get": self._get_split_ratio},
"status": {"get": self._get_status}
}
def _get_recuperator_total_UA(self):
return (
self.get_parameter("recuperator 1 UA")
+ self.get_parameter("recuperator 2 UA")
)
def _get_status(self):
return self.nw.status
def _calc_thermal_efficiency(self):
return (
self.nw.get_conn("e1").E.val_SI
/ self.nw.get_conn("h1").E.val_SI
)
def _set_split_ratio(self, value):
m = None
c10, c15 = self.nw.get_conn(["10", "15"])
if value is not None:
m=Ref(c15, value, 0)
c10.set_attr(m=m)
def _get_split_ratio(self):
c10, c15 = self.nw.get_conn(["10", "15"])
return c10.m.val / c15.m.val
def create_diagram(self, fluid_name):
"""Fluid property diagram with isolines suited for the cycle.
The default isoline range of the base class ends at 200 °C and
slightly above the critical pressure. The cycle however reaches
turbine inlet temperatures of 600 °C and pressures of more than
250 bar, so the temperature range and the pressure isolines are
extended (including the internal iteration bound p_max, which the
isotherms and isochores are calculated up to).
"""
diagram = FluidPropertyDiagram(fluid_name)
diagram.set_unit_system(self.nw.units)
diagram.set_isolines_subcritical(-20, 650)
diagram.set_isolines(
p=np.array([1, 5, 10, 25, 50, 100, 150, 200, 250, 300])
)
diagram.p_max = 300e5
diagram.calc_isolines()
return diagram
def solve_model(self, **kwargs):
self.solve_model_design(**kwargs)
def solve_design(self):
self.nw.solve(mode="design")
self.nw.assert_convergence()
def run_exergy_analysis(self):
from exerpy import ExergyAnalysis
E_F = {"inputs": ["h1"], "outputs": []}
E_P = {"inputs": ["e1"], "outputs": []}
self._ean = ExergyAnalysis.from_tespy(
self.nw, self.Tamb + 273.15, self.pamb * 1e5
)
# the cooler rejects heat to the ambient without a material cooling
# stream; it is a dissipative component (as in the publication), its
# exergy input is destroyed rather than delivered as a product
self._ean.components["Water cooler"].dissipative = True
self._ean.analyse(E_F=E_F, E_P=E_P)
return self._ean
Design simulation¶
Instantiating the class builds the network once, solve_model runs the design
simulation. The net power output is imposed at 100 MW.
model = SCO2Model()
model.solve_model()
indicators = [
"net power", "heat input", "thermal efficiency",
"turbine inlet temperature", "turbine inlet pressure",
"compressor inlet temperature", "compressor inlet pressure",
]
pd.Series(model.get_results(indicators), name="design point").round(4)
net power 100.0000
heat input 245.9743
thermal efficiency 0.4065
turbine inlet temperature 600.0000
turbine inlet pressure 250.0000
compressor inlet temperature 35.0000
compressor inlet pressure 75.0000
Name: design point, dtype: float64
Validation¶
The reference publication proves temperature, pressure as well as thermal and mechanical specific exergies for every connection. All deviations are close to zero on the state variables. The mechanical exergy has the largest relative deviation, which comes from a different reference state implementation rather than from the cycle model.
published = pd.DataFrame(
# label: T in °C, p in bar, e_T in kJ/kg, e_M in kJ/kg
[[ 1, 35.00, 75.00, 7.59, 198.46],
[ 2, 123.26, 258.40, 34.86, 218.09],
[ 3, 433.63, 257.00, 231.28, 217.95],
[ 4, 600.00, 250.00, 362.76, 217.23],
[ 5, 456.95, 77.95, 213.80, 198.80],
[ 6, 128.26, 75.15, 33.56, 198.48],
[ 10, 128.26, 75.15, 33.56, 198.48],
[ 11, 264.11, 257.51, 116.43, 218.00],
[ 12, 264.11, 257.51, 116.43, 218.00],
[ 13, 264.11, 257.51, 116.43, 218.00],
[ 14, 269.11, 76.94, 95.93, 198.68],
[ 15, 128.26, 75.15, 33.56, 198.48]],
columns=["label", "T", "p", "e_T", "e_M"],
).set_index("label")
tespy = pd.DataFrame({
"T": {int(c.label): c.T.val for c in model.nw.conns["object"]
if c.label.isdigit()},
"p": {int(c.label): c.p.val for c in model.nw.conns["object"]
if c.label.isdigit()},
}).sort_index()
comparison = published[["T", "p"]].join(
tespy, lsuffix=" published", rsuffix=" tespy"
)
comparison["ΔT in K"] = comparison["T tespy"] - comparison["T published"]
comparison["δp in %"] = (
(comparison["p tespy"] - comparison["p published"])
/ comparison["p published"] * 100
)
comparison.round(3)
| T published | p published | T tespy | p tespy | ΔT in K | δp in % | |
|---|---|---|---|---|---|---|
| label | ||||||
| 1 | 35.00 | 75.00 | 35.000 | 75.00 | 0.000 | 0.0 |
| 2 | 123.26 | 258.40 | 123.235 | 258.40 | -0.025 | 0.0 |
| 3 | 433.63 | 257.00 | 433.693 | 257.00 | 0.063 | 0.0 |
| 4 | 600.00 | 250.00 | 600.000 | 250.00 | 0.000 | 0.0 |
| 5 | 456.95 | 77.95 | 457.141 | 77.95 | 0.191 | 0.0 |
| 6 | 128.26 | 75.15 | 128.235 | 75.15 | -0.025 | 0.0 |
| 10 | 128.26 | 75.15 | 128.235 | 75.15 | -0.025 | 0.0 |
| 11 | 264.11 | 257.51 | 264.139 | 257.51 | 0.029 | 0.0 |
| 12 | 264.11 | 257.51 | 264.139 | 257.51 | 0.029 | 0.0 |
| 13 | 264.11 | 257.51 | 264.139 | 257.51 | 0.029 | 0.0 |
| 14 | 269.11 | 76.94 | 269.139 | 76.94 | 0.029 | 0.0 |
| 15 | 128.26 | 75.15 | 128.235 | 75.15 | -0.025 | 0.0 |
Ts diagram¶
We can make the cycle plot in the Ts diagram:
fig, ax = model.plot_Ts_diagram_matplotlib("01")
RuntimeWarning: divide by zero encountered in divide
RuntimeWarning: invalid value encountered in divide
RuntimeWarning: invalid value encountered in multiply
Sensitivity analysis¶
At the example of the paramters below a sensitivity analysis is carried out:
parameter |
design value |
range |
|---|---|---|
turbine inlet temperature |
600 °C |
450 - 700 °C |
compressor inlet temperature |
35 °C |
32 - 45 °C |
compressor inlet pressure |
75 bar |
70 - 75 bar |
# keep build time for docs low
n = 6 if os.getenv("TESPY_DOCS_BUILD") else 13
tit = model.sensitivity_analysis(
param_dict={"turbine inlet temperature": np.linspace(450, 700, n)},
result_param_list=["thermal efficiency", "heat input"],
)
cit = model.sensitivity_analysis(
param_dict={"compressor inlet temperature": np.linspace(32, 45, n)},
result_param_list=["thermal efficiency"],
)
cip = model.sensitivity_analysis(
param_dict={"compressor inlet pressure": np.linspace(70, 75, n)},
result_param_list=["thermal efficiency"],
)
# solve_model() re-solves with whatever parameters are currently set,
# so the design values are restated rather than assumed
model.solve_model(**{'turbine inlet temperature': 600, 'compressor inlet temperature': 35, 'compressor inlet pressure': 75})
The turbine inlet temperature has the stronges effect on efficiency. The compressor inlet conditions also have a quite high effect, mostly due to the fact, that the fluid near to the critical point, where the density is high and the compression work small. Moving away from that point reduces efficiency.
fig, ax = plt.subplots(1, 3, figsize=(12, 3.6))
ax[0].plot(tit["turbine inlet temperature"], tit["thermal efficiency"], "o-", color="tab:red")
ax[0].set_xlabel("turbine inlet temperature in °C")
ax[1].plot(cit["compressor inlet temperature"], cit["thermal efficiency"], "o-", color="tab:blue")
ax[1].set_xlabel("compressor inlet temperature in °C")
ax[1].axvline(31.0, ls="--", color="gray", lw=1)
ax[1].annotate(
"critical\ntemperature", (31.0, cit["thermal efficiency"].min()),
textcoords="offset points", xytext=(6, 4), fontsize=8, color="gray"
)
ax[2].plot(cip["compressor inlet pressure"], cip["thermal efficiency"], "o-", color="tab:green")
ax[2].set_xlabel("compressor inlet pressure in bar")
ax[2].axvline(73.8, ls="--", color="gray", lw=1)
ax[2].annotate(
"critical\npressure", (73.8, cip["thermal efficiency"].min()),
textcoords="offset points", xytext=(-38, 4), fontsize=8, color="gray"
)
for a in ax:
a.set_ylabel("thermal efficiency")
a.grid(alpha=.3)
plt.tight_layout()
plt.close(fig)
fig
Changing the split ratio¶
The original data from the reference publication assume a specific split ratio of the mass flow between the compressor and the recomplressor. This is done implicitly through the temperature reference of connection 11. We can run a sensitivity analysis over this ratio:
c10, c11, c15 = (model.nw.get_conn(x) for x in ("10", "11", "15"))
c11.set_attr(T=None) # release the matched merge temperature
result = model.sensitivity_analysis(
param_dict={"split ratio": np.linspace(0.15, 0.30, 6)},
result_param_list=["thermal efficiency", "T recompressed in", "T recuperated in", "status"],
)
result.set_index("split ratio").round(
{"thermal efficiency": 3, "status": 0,
"T recompressed in": 1, "T recuperated in": 1}
)
Invalid value for ttd_u: ttd_u = -9.769978344902142 below minimum value (0) at component Recuperator 1.
Invalid value for ttd_min: ttd_min = -9.769978344902142 below minimum value (0) at component Recuperator 1.
Invalid value for eff_cold: eff_cold = 1.0464114775961735 above maximum value (1) at component Recuperator 1.
Invalid value for eff_max: eff_max = 1.0464114775961735 above maximum value (1) at component Recuperator 1.
Invalid value for td_pinch: td_pinch = -9.769978344902142 below minimum value (0) at component Recuperator 1.
| thermal efficiency | T recompressed in | T recuperated in | status | |
|---|---|---|---|---|
| split ratio | ||||
| 0.15 | 0.392 | 264.1 | 161.7 | 0 |
| 0.18 | 0.395 | 264.1 | 176.9 | 0 |
| 0.21 | 0.399 | 264.1 | 196.9 | 0 |
| 0.24 | 0.402 | 264.1 | 223.9 | 0 |
| 0.27 | 0.406 | 264.1 | 261.1 | 0 |
| 0.30 | 0.410 | 264.1 | 312.8 | 1 |
We monitor the highest efficiency at the highest split ratio, but we also see
that the solution for that state is impossible. status==1 indicates a
vioalation of what is physically possible. In this case, the recuperator 1 has
a crossing of the two temperature profiles. We can prevent this, by imposing
the minimum pinch instead of the terminal temperature differences.
model.solve_model(**{"split ratio": 0.3})
pd.DataFrame({
label: {
"ttd upper in K": model.nw.get_comp(label).ttd_u.val,
"ttd lower in K": model.nw.get_comp(label).ttd_l.val,
"pinch in K": model.nw.get_comp(label).td_pinch.val,
"UA in MW/K": model.nw.get_comp(label).UA.val,
}
for label in ("Recuperator 1", "Recuperator 2")
}).round(3)
Invalid value for ttd_u: ttd_u = -9.769978344880201 below minimum value (0) at component Recuperator 1.
Invalid value for ttd_min: ttd_min = -9.769978344880201 below minimum value (0) at component Recuperator 1.
Invalid value for eff_cold: eff_cold = 1.0464114775960642 above maximum value (1) at component Recuperator 1.
Invalid value for eff_max: eff_max = 1.0464114775960642 above maximum value (1) at component Recuperator 1.
Invalid value for td_pinch: td_pinch = -9.769978344880201 below minimum value (0) at component Recuperator 1.
| Recuperator 1 | Recuperator 2 | |
|---|---|---|
| ttd upper in K | -9.77 | 18.537 |
| ttd lower in K | 5.00 | 5.000 |
| pinch in K | -9.77 | 5.000 |
| UA in MW/K | NaN | 19.044 |
We can also plot the temperature profiles:
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
model.plot_QT_diagram_matplotlib("Recuperator 1", ax[0])
model.plot_QT_diagram_matplotlib("Recuperator 2", ax[1])
for a in ax:
a.grid(alpha=.3)
plt.tight_layout()
plt.close(fig)
fig
So instead we switch the specifications to make use of td_pinch instead.
model.set_parameters(**{
"recuperator 1 ttd_l": None, "recuperator 2 ttd_l": None,
"recuperator 1 pinch": 5, "recuperator 2 pinch": 5
})
result = model.sensitivity_analysis(
param_dict={"split ratio": np.linspace(0.15, 0.30, 6)},
result_param_list=["thermal efficiency", "T recompressed in", "T recuperated in", "status"],
)
result.set_index("split ratio").round(
{"thermal efficiency": 3, "status": 0,
"T recompressed in": 1, "T recuperated in": 1}
)
| thermal efficiency | T recompressed in | T recuperated in | status | |
|---|---|---|---|---|
| split ratio | ||||
| 0.15 | 0.392 | 264.1 | 161.7 | 0 |
| 0.18 | 0.395 | 264.1 | 176.9 | 0 |
| 0.21 | 0.399 | 264.1 | 196.9 | 0 |
| 0.24 | 0.402 | 264.1 | 223.9 | 0 |
| 0.27 | 0.406 | 264.1 | 261.1 | 0 |
| 0.30 | 0.397 | 272.1 | 272.1 | 0 |
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
model.plot_QT_diagram_matplotlib("Recuperator 1", ax[0])
model.plot_QT_diagram_matplotlib("Recuperator 2", ax[1])
ax[0].set_title("Recuperator 1")
ax[1].set_title("Recuperator 2")
for a in ax:
a.grid(alpha=.3)
plt.tight_layout()
plt.close(fig)
fig
Recuperator size analysis¶
Instead of varying the split ratio, we carry out a sensitivity analysis on the sizing of the recuperators by imposing a lower or higher pinch value. For this we bring back the original split ratio defined by temperature equality and impose the lower terminal temperature difference. We look at recuperator 1 only at first:
model.set_parameters(
**{"split ratio": None, "recuperator 1 pinch": None, "recuperator 1 ttd_l": 5}
)
c11, c12 = model.nw.get_conn(["11", "12"])
c11.set_attr(T=Ref(c12, 1, 0))
model.solve_model()
n = 6 if os.getenv("TESPY_DOCS_BUILD") else 12
sizing = model.sensitivity_analysis(
param_dict={"recuperator 1 ttd_l": np.linspace(2, 18, n)},
result_param_list=["thermal efficiency", "recuperator 1 UA"],
)
fig, ax = plt.subplots(1, 2, figsize=(11, 3.8))
ax[0].plot(sizing["recuperator 1 ttd_l"], sizing["thermal efficiency"], "o-", color="tab:blue")
ax[0].set_xlabel("recuperator 1 terminal temperature difference in K")
ax[0].set_ylabel("thermal efficiency")
ax[1].plot(sizing["recuperator 1 UA"], sizing["thermal efficiency"], "o-", color="tab:red")
ax[1].set_xlabel("recuperator 1 UA in MW/K")
ax[1].set_ylabel("thermal efficiency")
for a in ax:
a.grid(alpha=.3)
plt.tight_layout()
plt.close(fig)
fig
We can see, that the efficiency gain is reduced with an increae in the recuperator’s heat transfer coefficient, meaning a low temeprature difference increases the necessary area of the heat exchanger by a lot.
Optimizing the two recuperators¶
After the sensitivity analysis we optimize both heat exchangers simultaneously. We want to achieve high efficiency and low overall UA at the same time. With the NSGA-II algorithm we can make a pareto optimization.
Note
The generation count is reduced during a documentation build. Run the notebook interactively for a complete output.
n_gen = 3 if os.getenv("TESPY_DOCS_BUILD") else 25
model.set_parameters(**{"recuperator 2 pinch": None})
log, result = model.optimize(
algorithm=NSGA2(pop_size=24),
termination=("n_gen", n_gen),
variables={
"recuperator 1 ttd_l": {"min": 2, "max": 18},
"recuperator 2 ttd_l": {"min": 2, "max": 18},
},
objective=["thermal efficiency", "recuperator UA"],
minimize_flags=[False, True],
kpi=["heat input"],
)
feasible = log.dropna(subset=["thermal efficiency", "recuperator UA"])
feasible.sort_values("thermal efficiency", ascending=False).head().round(3)
| recuperator 1 ttd_l | recuperator 2 ttd_l | thermal efficiency | recuperator UA | heat input | |
|---|---|---|---|---|---|
| 68 | 2.161 | 4.235 | 0.411 | 48.946 | 243.172 |
| 54 | 2.190 | 4.320 | 0.411 | 48.520 | 243.232 |
| 35 | 2.181 | 4.776 | 0.411 | 46.856 | 243.412 |
| 58 | 3.755 | 4.424 | 0.409 | 44.415 | 244.641 |
| 36 | 3.645 | 4.817 | 0.409 | 43.248 | 244.705 |
fig, ax = plt.subplots(figsize=(10.5, 5.8))
sc = ax.scatter(
feasible["recuperator UA"], feasible["thermal efficiency"],
c=feasible["recuperator 1 ttd_l"], cmap="viridis", s=26,
)
fig.colorbar(sc, label="recuperator 1 ttd in K")
ax.set_xlabel("recuperator combined UA in MW/K")
ax.set_ylabel("thermal efficiency")
ax.grid(alpha=.3)
plt.close(fig)
fig
