CHP plant & DHN design_Problem setting up the system #476

MikeSennis · 2024-01-17T19:30:06Z

MikeSennis
Jan 17, 2024

Hello,
I'm trying to model a CHP plant in TESPy, connected to a district heating network (DHN). The power plant consists of 2 gas turbines (of the same specs), the exhaust gases of which will heat DHN's water. Also, these exhaust gases will connect to a single stream through a merger ("Merge" component) before ending up in the DHN's heat exchanger (as shown in the uploaded figure). I've tried to build up my system in incremental steps, adding initially the first gas turbine, then the second one, and at the end, the flue gas merger and the DHN to my existing system. But for some reason, the following error message occurs when running the simulation:

Moreover, is there a way to add the 2nd gas turbine to the system without needing to set this gas turbine up from scratch (the two gas turbines are identical)?
Any help or suggestion with these issues would be highly appreciated!

System:

Code:

from tespy.networks import Network
from tespy.components import (
    DiabaticCombustionChamber, Turbine, Source, Sink, Compressor
)
from tespy.connections import Connection, Ref, Bus

nw = Network(p_unit="bar", T_unit="C", h_unit="kJ / kg")

# adding the 1st gas turbine to our existing network
cp1 = Compressor("compressor 1")
cc1 = DiabaticCombustionChamber("combustion chamber 1")
tu1 = Turbine("turbine 1")
air1 = Source("air source 1")
fuel1 = Source("fuel source 1")
fg1 = Sink("flue gas sink 1")

c8 = Connection(air1, "out1", cc1, "in1", label="8")
c9 = Connection(cc1, "out1", fg1, "in1", label="9")
c11 = Connection(fuel1, "out1", cc1, "in2", label="11")
nw.add_conns(c8, c9, c11)

cc1.set_attr(pr=0.98, eta=1, lamb=1.5, ti=None) # or ti=138.3899e6 W (LHV = 49.736 MJ/kg for natural gas)

c8.set_attr(
    p=0.93194, T=13.5,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions, m=122.185 [kg/s]
c11.set_attr(m=0.0001, p=Ref(c8, 1.05, 0), T=13.5, fluid={"CH4": 0.96, "ETHANE": 0.04, "H2": 0}) # fluid: natural gas -> composition by mass fractions

nw.solve(mode="design")

#cc.set_attr(ti=None)
c11.set_attr(m=None)
c11.set_attr(m=1)
nw.solve(mode="design")

c11.set_attr(m=None)
c11.set_attr(m=2.78249) 
nw.solve(mode="design")

cc1.set_attr(lamb=None)
c9.set_attr(T=1306.85)
nw.solve(mode="design")
nw.print_results()

nw.del_conns(c8, c9)
c7 = Connection(air1, "out1", cp1, "in1", label="7")
c8 = Connection(cp1, "out1", cc1, "in1", label="8")
c9 = Connection(cc1, "out1", tu1, "in1", label="9")
c10 = Connection(tu1, "out1", fg1, "in1", label="10")
nw.add_conns(c7, c8, c9, c10)

generator1 = Bus("generator 1")
generator1.add_comps(
    {"comp": tu1, "char": 1, "base": "component"},
    {"comp": cp1, "char": 1, "base": "bus"},
)
nw.add_busses(generator1)

cp1.set_attr(eta_s=0.9, pr=21.1)
tu1.set_attr(eta_s=0.89) 
c7.set_attr(
    m=None, p=0.93194, T=13.5,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions
c9.set_attr(m=124.967) 
c10.set_attr(p=Ref(c7, 1.03, 0))
c11.set_attr(p=None)
c11.set_attr(p=Ref(c8, 1.05, 0))
nw.solve("design")

c9.set_attr(m=None, T=1306.85) # or m=124.967
nw.solve("design")
nw.print_results()
el_eff1 = (abs(generator1.P.val) / cc1.ti.val) * 100 # Brayton cycle (electric) efficiency
print(f"Brayton cycle 1 (electric) efficiency = {round(el_eff1, 2)} %", "\n")
print(c10.fluid_data)

# adding the 2nd gas turbine to our existing network
cp2 = Compressor("compressor 2")
cc2 = DiabaticCombustionChamber("combustion chamber 2")
tu2 = Turbine("turbine 2")
air2 = Source("air source 2")
fuel2 = Source("fuel source 2")
fg2 = Sink("flue gas sink 2")

c13 = Connection(air2, "out1", cc2, "in1", label="13")
c14 = Connection(cc2, "out1", fg2, "in1", label="14")
c16 = Connection(fuel2, "out1", cc2, "in2", label="16")
nw.add_conns(c13, c14, c16)

cc2.set_attr(pr=0.98, eta=1, lamb=1.5, ti=None) # or ti=138.3899e6 W (LHV = 49.736 MJ/kg for natural gas)

c13.set_attr(
    p=0.93194, T=13.5,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions, m=122.185 [kg/s]
c16.set_attr(m=0.0001, p=Ref(c13, 1.05, 0), T=13.5, fluid={"CH4": 0.96, "ETHANE": 0.04, "H2": 0}) # fluid: natural gas -> composition by mass fractions

nw.solve(mode="design")

#cc.set_attr(ti=None)
c16.set_attr(m=None)
c16.set_attr(m=1) 
nw.solve(mode="design")

c16.set_attr(m=None)
c16.set_attr(m=2.78249) 
nw.solve(mode="design")

cc2.set_attr(lamb=None)
c14.set_attr(T=1306.85)
nw.solve(mode="design")
nw.print_results()

nw.del_conns(c13, c14)
c12 = Connection(air2, "out1", cp2, "in1", label="12")
c13 = Connection(cp2, "out1", cc2, "in1", label="13")
c14 = Connection(cc2, "out1", tu2, "in1", label="14")
c15 = Connection(tu2, "out1", fg2, "in1", label="15")
nw.add_conns(c12, c13, c14, c15)

generator2 = Bus("generator 2")
generator2.add_comps(
    {"comp": tu2, "char": 1, "base": "component"},
    {"comp": cp2, "char": 1, "base": "bus"},
)
nw.add_busses(generator2)

cp2.set_attr(eta_s=0.9, pr=21.1)
tu2.set_attr(eta_s=0.89) 
c12.set_attr(
    m=None, p=0.93194, T=13.5,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions, m=122.185 
c14.set_attr(m=124.967)
c15.set_attr(p=Ref(c12, 1.03, 0))
c16.set_attr(p=None)
c16.set_attr(p=Ref(c13, 1.05, 0))
nw.solve("design")

c14.set_attr(m=None, T=1306.85) # or m=124.967
nw.solve("design")
nw.print_results()
el_eff2 = (abs(generator2.P.val) / cc2.ti.val) * 100 # Brayton cycle (electric) efficiency
print(f"Brayton cycle 2 (electric) efficiency = {round(el_eff2, 2)} %", "\n")
print(c15.fluid_data)

# adding the flue gas merger and the district heating network to our existing network
from tespy.components import (
    Merge, CycleCloser, Pipe, Pump, Valve, SimpleHeatExchanger, HeatExchanger
)

# central heating plant (CHP plant)
me = Merge('flue gas merger') # Power plant's mixing chamber
hs = HeatExchanger('heat source')  # Power plant's heat exchanger
cc = CycleCloser('cycle closer')
fp = Pump('feed pump')
fg = Sink('flue gas sink')

# consumer
cons = SimpleHeatExchanger('consumer')
val = Valve('control valve')

# pipes
pipe_feed = Pipe('feed pipe')
pipe_return = Pipe('return pipe')

nw.del_conns(c10, c15)
# flue gas loop
c10 = Connection(tu1, "out1", me, "in1", label="10")
c15 = Connection(tu2, "out1", me, "in2", label="15")
c17 = Connection(me, "out1", hs, 'in1', label="17") # in1 for heat exchanger's hot side inlet
c18 = Connection(hs, 'out1', fg, 'in1', label='18') # out1 for heat exchanger's hot side outlet
nw.add_conns(c10, c15, c17, c18)

# main (closed) cycle: water cycle
c0 = Connection(cc, "out1", hs, "in2", label="0")
c1 = Connection(hs, "out2", fp, "in1", label="1")
c2 = Connection(fp, "out1", pipe_feed, "in1", label="2")
c3 = Connection(pipe_feed, "out1", cons, "in1", label="3")
c4 = Connection(cons, "out1", val, "in1", label="4")
c5 = Connection(val, "out1", pipe_return, "in1", label="5")
c6 = Connection(pipe_return, "out1", cc, "in1", label="6")
nw.add_conns(c0, c1, c2, c3, c4, c5, c6)

# parameters' definition
hs.set_attr(pr2=1)
cons.set_attr(pr=0.98)
fp.set_attr(eta_s=0.75)
pipe_feed.set_attr(Q=-(0.03*128.775e+06), L=17140, D=0.6, pr=0.98)  # assuming 4.8 % total thermal losses in the pipes (128.775 MWth from analysis in Excel)
pipe_return.set_attr(Q=-(0.018*128.775e+06), L=32216, D=2*0.45, pr=0.98)  # assuming 4.8 % total thermal losses in the pipes (128.775 MWth from analysis in Excel)

c10.set_attr(m=126.6, h=1053, p=0.9599, fluid={"H2O": 0.04895, "CO2": 0.06102, "N2": 0.73841, "O2": 0.139, "Ar": 0.01262})
                                        # fluid: flue gas -> composition by mass fractions
c15.set_attr(m=126.6, h=1053, p=0.9599, fluid={"H2O": 0.04895, "CO2": 0.06102, "N2": 0.73841, "O2": 0.139, "Ar": 0.01262})
                                        # fluid: flue gas -> composition by mass fractions
c18.set_attr(T=120, p=0.93194) 
c1.set_attr(T=115, p=23, fluid={'HEOS::Water': 1})
c2.set_attr(p=25)
c0.set_attr(T=65)

nw.solve(mode="design")
nw.print_results()
print(f"DHN mass flow rate = {round(c0.m.val, 3)} kg/s")
print(f"Thermal power given from CHP plant to (primary) DHN = {round(abs(hs.Q.val)/1e6, 3)} MWth")
print(f"Thermal power given from primary DHN to secondary DHN (consumers) = {round(abs(cons.Q.val)/1e6, 3)} MWth")
print(f"Thermal efficiency of (primary) DHN = {round((abs(cons.Q.val)/abs(hs.Q.val)) * 100, 2)} %")

Answered by fwitte

Jan 21, 2024

And for the second option, in case you have the identical gas turbine setup:

from tespy.networks import Network
from tespy.components import (
    DiabaticCombustionChamber, Turbine, Source, Sink, Compressor, Subsystem
)
from tespy.connections import Connection, Ref, Bus


class GasTurbineSubsystem(Subsystem):
    """Class documentation"""

    def create_comps(self):
        """Create the subsystem's components."""
        self.comps['compressor'] = Compressor(f'{self.label}-compressor')
        self.comps['combustion'] = DiabaticCombustionChamber(f'{self.label}-combustion chamber')
        self.comps['gas'] = Source(f'{self.label}-gas intake')
        self.comps['air'] = Source(f'{self.l…

View full answer

fwitte · 2024-01-21T16:55:55Z

fwitte
Jan 21, 2024
Maintainer

Hi @MikeSennis and thanks for reaching out! The error showing for you is actually a bug, I will fix it asap.

Can you give me a couple more insights, if the two gas turbines always will be operated identically? In that case you could also remove the second one from the system and just double the flue gas mass flow...

If both gas turbines should not be identical, you could use the subsystem class of tespy to build a gas turbine as a subsystem and include that as desired, or you go even beyond and build your own class to handle that. See below for an example

from tespy.networks import Network
from tespy.components import (
    DiabaticCombustionChamber, Turbine, Source, Sink, Compressor, Subsystem
)
from tespy.connections import Connection, Ref, Bus


class GasTurbineSubsystem(Subsystem):
    """Class documentation"""

    def create_comps(self):
        """Create the subsystem's components."""
        self.comps['compressor'] = Compressor(f'{self.label}-compressor')
        self.comps['combustion'] = DiabaticCombustionChamber(f'{self.label}-combustion chamber')
        self.comps['gas'] = Source(f'{self.label}-gas intake')
        self.comps['air'] = Source(f'{self.label}-air intake')
        self.comps['turbine'] = Turbine(f'{self.label}-turbine')

    def create_conns(self):
        """Define the subsystem's connections."""
        self.conns['c1'] = Connection(
            self.comps['air'], 'out1', self.comps['compressor'], 'in1',
            label=f"{self.label}-1"
        )
        self.conns['c2'] = Connection(
            self.comps['compressor'], 'out1', self.comps['combustion'], 'in1',
            label=f"{self.label}-2"
        )
        self.conns['c3'] = Connection(
            self.comps['combustion'], 'out1', self.comps['turbine'], 'in1',
            label=f"{self.label}-3"
        )
        self.conns['c4'] = Connection(
            self.comps['gas'], 'out1', self.comps['combustion'], 'in2',
            label=f"{self.label}-4",
            # you can also include specifications here
            # if you want to modify them on creating the class instance, then
            # you have to modify the __init__ method and set up respective
            # logic and variables
            p=Ref(self.conns["c2"], 1.05, 0)
        )

nw = Network(p_unit="bar", T_unit="C", h_unit="kJ / kg")

# adding the 1st gas turbine to our existing network
gt1 = GasTurbineSubsystem("gt1")
fg1 = Sink("flue gas sink 1")

c10 = Connection(gt1.comps["turbine"], "out1", fg1, "in1", label="10")

nw.add_subsys(gt1)
nw.add_conns(c10)

generator1 = Bus("generator 1")
generator1.add_comps(
    {"comp": gt1.comps["turbine"], "char": 1, "base": "component"},
    {"comp": gt1.comps["compressor"], "char": 1, "base": "bus"},
)
nw.add_busses(generator1)

gt1.comps["combustion"].set_attr(pr=0.98, eta=1, lamb=1.5, ti=None) # or ti=138.3899e6 W (LHV = 49.736 MJ/kg for natural gas)

gt1.conns["c1"].set_attr(
    p=0.93194, T=13.5,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions, m=122.185 [kg/s]
gt1.conns["c4"].set_attr(m=1, T=13.5, fluid={"CH4": 0.96, "ETHANE": 0.04, "H2": 0}) # fluid: natural gas -> composition by mass fractions
c10.set_attr(p=Ref(gt1.conns["c1"], 1.03, 0))

gt1.comps["compressor"].set_attr(eta_s=0.9, pr=21.1)
gt1.comps["turbine"].set_attr(eta_s=0.89)

nw.solve(mode="design")

# adding the 2nd gas turbine to our existing network
gt2 = GasTurbineSubsystem("gt2")
fg2 = Sink("flue gas sink 2")

c20 = Connection(gt2.comps["turbine"], "out1", fg2, "in1", label="20")

nw.add_subsys(gt2)
nw.add_conns(c20)

generator2 = Bus("generator 2")
generator2.add_comps(
    {"comp": gt2.comps["turbine"], "char": 1, "base": "component"},
    {"comp": gt2.comps["compressor"], "char": 1, "base": "bus"},
)
nw.add_busses(generator2)

gt2.comps["combustion"].set_attr(pr=0.98, eta=1, lamb=1.5, ti=None) # or ti=138.3899e6 W (LHV = 49.736 MJ/kg for natural gas)

gt2.conns["c1"].set_attr(
    p=0.93194, T=13.5,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions, m=122.185 [kg/s]
gt2.conns["c4"].set_attr(m=1, T=13.5, fluid={"CH4": 0.96, "ETHANE": 0.04, "H2": 0}) # fluid: natural gas -> composition by mass fractions
c20.set_attr(p=Ref(gt2.conns["c1"], 1.03, 0))

gt2.comps["compressor"].set_attr(eta_s=0.9, pr=21.1)
gt2.comps["turbine"].set_attr(eta_s=0.89)

nw.solve(mode="design")
nw.print_results()

# adding the flue gas merger and the district heating network to our existing network
from tespy.components import (
    Merge, CycleCloser, Pipe, Pump, Valve, SimpleHeatExchanger, HeatExchanger
)

# central heating plant (CHP plant)
me = Merge('flue gas merger') # Power plant's mixing chamber
hs = SimpleHeatExchanger('heat source')  # Power plant's heat exchanger
fg = Sink('flue gas sink')

nw.del_conns(c10, c20)

c10 = Connection(gt1.comps["turbine"], "out1", me, "in1", label="10")
c20 = Connection(gt2.comps["turbine"], "out1", me, "in2", label="20")
c17 = Connection(me, "out1", hs, "in1", label="17")
c18 = Connection(hs, "out1", fg, "in1", label="18")

nw.add_conns(c10, c20, c17, c18)

c18.set_attr(p=Ref(gt1.conns["c1"], 1, 0), T=120)
hs.set_attr(pr=.99)

gt1.conns["c4"].set_attr(m=None)
gt2.conns["c4"].set_attr(m=None)

gt1.conns["c3"].set_attr(m=125)
gt2.conns["c3"].set_attr(m=125)

nw.solve("design")

nw.del_conns(c17, c18)

hs = HeatExchanger("waste heat recovery")

cc = CycleCloser('cycle closer')
fp = Pump('feed pump')
fg = Sink('flue gas sink')

# consumer
cons = SimpleHeatExchanger('consumer')
val = Valve('control valve')

# pipes
pipe_feed = Pipe('feed pipe')
pipe_return = Pipe('return pipe')

# flue gas loop
c17 = Connection(me, "out1", hs, 'in1', label="17") # in1 for heat exchanger's hot side inlet
c18 = Connection(hs, 'out1', fg, 'in1', label='18') # out1 for heat exchanger's hot side outlet
nw.add_conns(c17, c18)

# main (closed) cycle: water cycle
c0 = Connection(cc, "out1", hs, "in2", label="0")
c1 = Connection(hs, "out2", fp, "in1", label="1")
c2 = Connection(fp, "out1", pipe_feed, "in1", label="2")
c3 = Connection(pipe_feed, "out1", cons, "in1", label="3")
c4 = Connection(cons, "out1", val, "in1", label="4")
c5 = Connection(val, "out1", pipe_return, "in1", label="5")
c6 = Connection(pipe_return, "out1", cc, "in1", label="6")
nw.add_conns(c0, c1, c2, c3, c4, c5, c6)

# parameters' definition
hs.set_attr(pr2=1, pr1=1)
cons.set_attr(pr=0.98)
fp.set_attr(eta_s=0.75)
pipe_feed.set_attr(Q=-(0.03*128.775e+06), L=17140, D=0.6, pr=0.98)  # assuming 4.8 % total thermal losses in the pipes (128.775 MWth from analysis in Excel)
pipe_return.set_attr(Q=-(0.018*128.775e+06), L=32216, D=2*0.45, pr=0.98)  # assuming 4.8 % total thermal losses in the pipes (128.775 MWth from analysis in Excel)

# c10.set_attr(m=126.6, h=1053, p=0.9599, fluid={"H2O": 0.04895, "CO2": 0.06102, "N2": 0.73841, "O2": 0.139, "Ar": 0.01262})
#                                         # fluid: flue gas -> composition by mass fractions
# c15.set_attr(m=126.6, h=1053, p=0.9599, fluid={"H2O": 0.04895, "CO2": 0.06102, "N2": 0.73841, "O2": 0.139, "Ar": 0.01262})
#                                         # fluid: flue gas -> composition by mass fractions
c18.set_attr(T=120, p=0.93194)
# using incompressible water for district heating should be faster
c1.set_attr(T=115, p=23, fluid={'INCOMP::Water': 1})
c2.set_attr(p=25)
c0.set_attr(T=65)

# set starting values for incompressible fluid connections (starting value generator is not working well with INCOMPs)
[c.set_attr(h0=100) for c in [c0, c1, c2, c3, c4, c5, c6]]

nw.solve(mode="design")
nw.print_results()
print(f"DHN mass flow rate = {round(c0.m.val, 3)} kg/s")
print(f"Thermal power given from CHP plant to (primary) DHN = {round(abs(hs.Q.val)/1e6, 3)} MWth")
print(f"Thermal power given from primary DHN to secondary DHN (consumers) = {round(abs(cons.Q.val)/1e6, 3)} MWth")
print(f"Thermal efficiency of (primary) DHN = {round((abs(cons.Q.val)/abs(hs.Q.val)) * 100, 2)} %")

2 replies

fwitte Jan 21, 2024
Maintainer

The fix for the bug is applied here: #481

fwitte Jan 21, 2024
Maintainer

And for the second option, in case you have the identical gas turbine setup:

from tespy.networks import Network
from tespy.components import (
    DiabaticCombustionChamber, Turbine, Source, Sink, Compressor, Subsystem
)
from tespy.connections import Connection, Ref, Bus


class GasTurbineSubsystem(Subsystem):
    """Class documentation"""

    def create_comps(self):
        """Create the subsystem's components."""
        self.comps['compressor'] = Compressor(f'{self.label}-compressor')
        self.comps['combustion'] = DiabaticCombustionChamber(f'{self.label}-combustion chamber')
        self.comps['gas'] = Source(f'{self.label}-gas intake')
        self.comps['air'] = Source(f'{self.label}-air intake')
        self.comps['turbine'] = Turbine(f'{self.label}-turbine')

    def create_conns(self):
        """Define the subsystem's connections."""
        self.conns['c1'] = Connection(
            self.comps['air'], 'out1', self.comps['compressor'], 'in1',
            label=f"{self.label}-1"
        )
        self.conns['c2'] = Connection(
            self.comps['compressor'], 'out1', self.comps['combustion'], 'in1',
            label=f"{self.label}-2"
        )
        self.conns['c3'] = Connection(
            self.comps['combustion'], 'out1', self.comps['turbine'], 'in1',
            label=f"{self.label}-3"
        )
        self.conns['c4'] = Connection(
            self.comps['gas'], 'out1', self.comps['combustion'], 'in2',
            label=f"{self.label}-4",
            # you can also include specifications here
            # if you want to modify them on creating the class instance, then
            # you have to modify the __init__ method and set up respective
            # logic and variables
            p=Ref(self.conns["c2"], 1.05, 0)
        )

nw = Network(p_unit="bar", T_unit="C", h_unit="kJ / kg")

# adding the 1st gas turbine to our existing network
gt1 = GasTurbineSubsystem("gt1")
fg1 = Sink("flue gas sink 1")

c10 = Connection(gt1.comps["turbine"], "out1", fg1, "in1", label="10")

nw.add_subsys(gt1)
nw.add_conns(c10)

generator1 = Bus("generator 1")
generator1.add_comps(
    {"comp": gt1.comps["turbine"], "char": 1, "base": "component"},
    {"comp": gt1.comps["compressor"], "char": 1, "base": "bus"},
)
nw.add_busses(generator1)

gt1.comps["combustion"].set_attr(pr=0.98, eta=1, lamb=1.5, ti=None) # or ti=138.3899e6 W (LHV = 49.736 MJ/kg for natural gas)

gt1.conns["c1"].set_attr(
    p=0.93194, T=13.5,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions, m=122.185 [kg/s]
gt1.conns["c4"].set_attr(m=1, T=13.5, fluid={"CH4": 0.96, "ETHANE": 0.04, "H2": 0}) # fluid: natural gas -> composition by mass fractions
c10.set_attr(p=Ref(gt1.conns["c1"], 1.03, 0))

gt1.comps["compressor"].set_attr(eta_s=0.9, pr=21.1)
gt1.comps["turbine"].set_attr(eta_s=0.89)

nw.solve(mode="design")

nw.print_results()

# adding the flue gas merger and the district heating network to our existing network
from tespy.components import (
    Merge, CycleCloser, Pipe, Pump, Valve, SimpleHeatExchanger, HeatExchanger
)

# central heating plant (CHP plant)
# instead of a merge we can put a source with the same fluid properties but double mass flow
me = Source('flue gas merger') # Power plant's mixing chamber
hs = SimpleHeatExchanger('heat source')  # Power plant's heat exchanger
fg = Sink('flue gas sink')

c17 = Connection(me, "out1", hs, "in1", label="17")
c18 = Connection(hs, "out1", fg, "in1", label="18")

nw.add_conns(c17, c18)

c10.set_attr(p=None)
c17.set_attr(p=Ref(c10, 1, 0), h=Ref(c10, 1, 0), m=Ref(c10, 2, 0), fluid={fluid: x for fluid, x in c10.fluid.val.items()})
c18.set_attr(p=Ref(gt1.conns["c1"], 1, 0), T=120)
hs.set_attr(pr=.99)


gt1.conns["c4"].set_attr(m=None)
gt1.conns["c3"].set_attr(m=125)

nw.solve("design")

# instead of directly specifying the fluid composition on conection 17
# (it may change due to specifications inside the gas turbine), this function
# will ensure idetical fluid composition on connection 17 and 10.

from tespy.tools.helpers import UserDefinedEquation

def fluid_ude(self):
    c1, c2 = self.conns
    fluid = self.params["fluid"]
    return c1.fluid.val[fluid] - c2.fluid.val[fluid]

def fluid_ude_jacobian(self):
    c1, c2 = self.conns
    fluid = self.params["fluid"]
    if fluid in c1.fluid.is_var:
        self.jacobian[c1.fluid.J_col[fluid]] = 1
    if fluid in c2.fluid.is_var:
        self.jacobian[c2.fluid.J_col[fluid]] = -1

for fluid in c10.fluid.val:
    f_ude = UserDefinedEquation(fluid, fluid_ude, fluid_ude_jacobian, [c10, c17], {"fluid": fluid})
    nw.add_ude(f_ude)

c17.set_attr(fluid={fluid: None for fluid in c17.fluid.is_set})
nw.solve("design")

Answer selected by fwitte

MikeSennis · 2024-01-24T11:13:28Z

MikeSennis
Jan 24, 2024
Author

Hello @fwitte and thank you for your time and useful recommendations, they helped a lot!
My intention is for the two gas turbines to operate identically (so removing the second one and doubling the flue gas mass flow does the job) but of course, setting up the system to satisfy the more generic case (two not identical gas turbines) gives more flexibility and generally makes the code smarter.
I have two further questions. Firstly, I'm not sure I understand the difference between setting directly the fluid composition of connection 17 identical to that of connection 10 (c17.set_attr(fluid={fluid: x for fluid, x in c10.fluid.val.items()})) and on the other hand, using the "UserDefinedEquation" class to ensure identical fluid compositions for these connections. Both methods seem to give always the desired result.
Secondly, when using the "DiabaticCombustionChamber" class always the following error message occurs:

Is it because I always set the combustion thermal efficiency ("eta") equal to 1?
Thanks again for your time.

3 replies

fwitte Jan 24, 2024
Maintainer

Firstly, I'm not sure I understand the difference between setting directly the fluid composition of connection 17 identical to that of connection 10 (c17.set_attr(fluid={fluid: x for fluid, x in c10.fluid.val.items()})) and on the other hand, using the "UserDefinedEquation" class to ensure identical fluid compositions for these connections. Both methods seem to give always the desired result.

In case you run a simulation like this, the fluid compositions will be different.

With the first setup and the given value for lamb, you get the same values

nw.solve("design")

print(c10.fluid.val)
print(c17.fluid.val)

If you now change the lamb value of the gt, you do not know the fluid composition of the flue gas prior simulation, thus you have different values for fluid composition on c10 and c17, see:

gt1.comps["combustion"].set_attr(lamb=2) # you do not know the resulting fluid composition

nw.solve("design")

print(c10.fluid.val)
print(c17.fluid.val)

Using the user defined equation values are identical:

c17.set_attr(fluid={fluid: None for fluid in c17.fluid.is_set})
nw.solve("design")
print(c10.fluid.val)
print(c17.fluid.val)

Secondly, when using the "DiabaticCombustionChamber" class always the following error message occurs: [image]

That is a numerical issue with rounding errors triggering the warning. You can disable the warning, by reducing the threshold for maximum value of Q_loss like so:

self.comps['combustion'].Q_loss.set_attr(max_val=1e-3)

Best!

MikeSennis Jan 25, 2024
Author

Hello @fwitte, I think I get it now. Thanks for the explanation!
Best regards!

MikeSennis Feb 5, 2024
Author

Hello @fwitte,
I have a follow-up question regarding the above network configuration (the two gas turbines CHP plant connected to a DHN). I'm trying to model the off-design system performance and in the ideal case, I want the turbine inlet temperature (TIT), the air mass flow rate (ma), and the fuel mass flow rate (mf) to be variable (simultaneously). But I can't find a way for even one of those parameters to be set as a variable in the off-design case because I think that the available off-design parameters for the "DiabaticCombustionChamber" component are too few. Thus, I want to ask if there is a way to overcome this obstacle so that some of these parameters (even one or two of the three) can become variable in the off-design system simulation.
Thanks again for your time.

Code:

# For an energy system with two or more identical gas turbines
from tespy.networks import Network
from tespy.components import (
    DiabaticCombustionChamber, Turbine, Source, Sink, Compressor, Subsystem
)
from tespy.connections import Connection, Ref, Bus


class GasTurbineSubsystem(Subsystem):
    """Class documentation"""

    def create_comps(self):
        """Create the subsystem's components."""
        self.comps['compressor'] = Compressor(f'{self.label}-compressor')
        self.comps['combustion chamber'] = DiabaticCombustionChamber(f'{self.label}-combustion chamber')
        # self.comps['combustion chamber'].Q_loss.set_attr(max_val=1e-3)
        self.comps['fuel'] = Source(f'{self.label}-fuel intake')
        self.comps['air'] = Source(f'{self.label}-air intake')
        self.comps['turbine'] = Turbine(f'{self.label}-turbine')

    def create_conns(self):
        """Define the subsystem's connections."""
        self.conns['c1'] = Connection(
            self.comps['air'], 'out1', self.comps['compressor'], 'in1',
            label=f"{self.label}-1"
        )
        self.conns['c2'] = Connection(
            self.comps['compressor'], 'out1', self.comps['combustion chamber'], 'in1',
            label=f"{self.label}-2"
        )
        self.conns['c3'] = Connection(
            self.comps['combustion chamber'], 'out1', self.comps['turbine'], 'in1',
            label=f"{self.label}-3"
        )
        self.conns['c4'] = Connection(
            self.comps['fuel'], 'out1', self.comps['combustion chamber'], 'in2',
            label=f"{self.label}-4",
            # you can also include specifications here
            # if you want to modify them on creating the class instance, then
            # you have to modify the __init__ method and set up respective
            # logic and variables
            p=Ref(self.conns['c2'], 1.05, 0),
            T=Ref(self.conns['c1'], 1, 0)
        )

nw = Network(p_unit="bar", T_unit="C", h_unit="kJ / kg")

# ambient conditions
p_ambient = 1.01325
T_ambient = 15

# adding the 1st gas turbine to our existing network
gt1 = GasTurbineSubsystem("gt1")
fg1 = Sink("flue gas sink 1")

c10 = Connection(gt1.comps["turbine"], "out1", fg1, "in1", label="10")

nw.add_subsys(gt1)
nw.add_conns(c10)

generator1 = Bus("generator 1")
generator1.add_comps(
    {"comp": gt1.comps["turbine"], "char": 1, "base": "component"},
    {"comp": gt1.comps["compressor"], "char": 1, "base": "bus"},
)
nw.add_busses(generator1)

gt1.comps["combustion chamber"].set_attr(pr=0.98, eta=1, lamb=1.5, ti=None)
                                # or ti=149.5323e6 W (LHV = 49.7365 MJ/kg for natural gas)

gt1.conns["c1"].set_attr(
    p=p_ambient, T=T_ambient,
    fluid={"N2": 0.755, "O2": 0.2315, "Ar": 0.0129, "CO2": 0.0006}
) # fluid: (dry) air -> composition by mass fractions, m=132.496 [kg/s], ISO conditions (design point)
gt1.conns["c4"].set_attr(m=3.00649, fluid={"CH4": 0.96, "ETHANE": 0.04, "H2": 0})
                                        # fluid: natural gas -> composition by mass fractions
c10.set_attr(p=Ref(gt1.conns["c1"], 1.03, 0))

gt1.comps["compressor"].set_attr(eta_s=0.9, pr=21.1)
gt1.comps["turbine"].set_attr(eta_s=0.89)

nw.solve(mode="design")

gt1.comps["combustion chamber"].set_attr(lamb=None)
gt1.conns["c3"].set_attr(m=135.502) 
nw.solve("design")

gt1.conns["c3"].set_attr(m=None, T=1306.85) # or m=135.502 
nw.solve(mode="design")

# adding the flue gas merger and the district heating network to our existing network
from tespy.components import (
    CycleCloser, Pipe, Pump, Valve, SimpleHeatExchanger, HeatExchanger
)

# central heating plant (CHP plant)
# instead of a merge we can put a source with the same fluid properties but double (gas turbine's 1) flue gas mass flow
# this model (method) works only for identical gas turbines (with the same specs)
me = Source('flue gas merger') # Power plant's mixing chamber
hs = HeatExchanger('waste heat recovery') # Power plant's heat exchanger
cc = CycleCloser('cycle closer')
fp = Pump('feed pump')
fg = Sink('flue gas sink')

# consumer
cons = SimpleHeatExchanger('consumer')
val = Valve('control valve')

# pipes
pipe_feed = Pipe('feed pipe')
pipe_return = Pipe('return pipe')

# flue gas loop
c17 = Connection(me, "out1", hs, "in1", label="17") # in1 for heat exchanger's hot side inlet
c18 = Connection(hs, "out1", fg, "in1", label="18") # out1 for heat exchanger's hot side outlet
nw.add_conns(c17, c18)

# main (closed) cycle: water cycle
c0 = Connection(cc, "out1", hs, "in2", label="0")
c1 = Connection(hs, "out2", fp, "in1", label="1")
c2 = Connection(fp, "out1", pipe_feed, "in1", label="2")
c3 = Connection(pipe_feed, "out1", cons, "in1", label="3")
c4 = Connection(cons, "out1", val, "in1", label="4")
c5 = Connection(val, "out1", pipe_return, "in1", label="5")
c6 = Connection(pipe_return, "out1", cc, "in1", label="6")
nw.add_conns(c0, c1, c2, c3, c4, c5, c6)

# parameters' definition
hs.set_attr(pr2=1)
cons.set_attr(pr=0.98)
fp.set_attr(eta_s=0.75)
pipe_feed.set_attr(Q=-(0.03*145.526e+06), pr=0.98, L=17140, D=0.6, ks=None)
            # assuming 4.8 % total thermal losses in the pipes (145.526 MWth from the first run),
            # ks=0.000045: pipe's (absolute) roughness in meters for Commercial or Welded Steel (or assuming pr=0.98)
pipe_return.set_attr(Q=-(0.018*145.526e+06), pr=0.98, L=32216, D=2*0.45, ks=None)
            # assuming 4.8 % total thermal losses in the pipes (145.526 MWth from the first run), pipe's roughness in meters
            # ks=0.000045: pipe's (absolute) roughness in meters for Commercial or Welded Steel (or assuming pr=0.98)

c17.set_attr(p=Ref(c10, 1, 0), h=Ref(c10, 1, 0), m=Ref(c10, 2, 0), fluid={fluid: x for fluid, x in c10.fluid.val.items()})
                                                                            # x (float): Gas component mass fraction specification
c18.set_attr(p=Ref(gt1.conns["c1"], 1, 0), T=120) # assuming that exhaust gases' exit total pressure is equal to ambient (static)
                                                    # pr1(hs) ~= 0,97
                                                    # we could also assume pr1(hs) ~= 1 (pexit/pamb = 1,03 as in GasTurb)
c1.set_attr(T=115, p=23, fluid={'HEOS::Water': 1}) # HEOS backend: uses IAPWS-95 formulation (very accurate), which is slower,
                            # but more accurate than the IAPWS-IF97 formulation (slightly more accurate)
                            # using incompressible water for district heating should be faster (but very inaccurate)
c2.set_attr(p=25)
c0.set_attr(T=65)

# setting starting values for incompressible fluid connections (starting value generator is not working well with INCOMP)
# [c.set_attr(h0=100) for c in [c0, c1, c2, c3, c4, c5, c6]]

nw.solve("design")
print(c10.fluid.val)
print(c17.fluid.val)

# instead of directly specifying the fluid composition on connection 17
# (it may change due to specifications inside the gas turbine), this function
# will ensure identical fluid composition on connections 17 and 10.

from tespy.tools.helpers import UserDefinedEquation

def fluid_ude(self): # function holding the equation to be applied: equation definition
    c1, c2 = self.conns
    fluid = self.params["fluid"]
    return c1.fluid.val[fluid] - c2.fluid.val[fluid]
        # returns the residual value of the equation

def fluid_ude_jacobian(self): # function holding the calculation of the Jacobian
    c1, c2 = self.conns
    fluid = self.params["fluid"]
    if fluid in c1.fluid.is_var:
        self.jacobian[c1.fluid.J_col[fluid]] = 1
    if fluid in c2.fluid.is_var:
        self.jacobian[c2.fluid.J_col[fluid]] = -1

for fluid in c10.fluid.val:
    f_ude = UserDefinedEquation(fluid, fluid_ude, fluid_ude_jacobian, [c10, c17], {"fluid": fluid})
    nw.add_ude(f_ude)

c17.set_attr(fluid={fluid: None for fluid in c17.fluid.is_set})

nw.solve("design")
nw.save("CHP & DHN_identical GTs_design point")
nw.print_results()
print(f"Brayton cycle 1 net electric power output = {round(abs(generator1.P.val)/1e6, 3)} MWe", "\n")
el_eff1 = (abs(generator1.P.val) / gt1.comps["combustion chamber"].ti.val) * 100 # Brayton cycle (electric) efficiency
print(f"Brayton cycle 1 (electric) efficiency = {round(el_eff1, 2)} %", "\n")
print(c10.fluid.val)
print(c17.fluid.val)
print(f"DHN mass flow rate = {round(c0.m.val, 3)} kg/s")
print(f"Thermal power given from CHP plant to (primary) DHN = {round(abs(hs.Q.val)/1e6, 3)} MWth")
print(f"Thermal power given from primary DHN to secondary DHN (consumers) = {round(abs(cons.Q.val)/1e6, 3)} MWth")
print(f"Thermal efficiency of (primary) DHN = {round((abs(cons.Q.val)/abs(hs.Q.val)) * 100, 2)} %", "\n")

# off-design performance
gt1.comps["compressor"].set_attr(design=["pr", "eta_s"], offdesign=["char_map_pr", "char_map_eta_s"])
gt1.comps["turbine"].set_attr(design=["eta_s"], offdesign=["eta_s_char", "cone"]) # Stodola's cone law?
fp.set_attr(design=["eta_s"], offdesign=["eta_s_char"])
#gt1.conns["c4"].set_attr(design=["m"]) # because I want variable fuel mass flow rate
#gt1.conns["c2"].set_attr(offdesign=["m"])
gt1.comps["combustion chamber"].set_attr(offdesign=["lamb"]) # "ti" or "lamb"
#gt1.conns["c3"].set_attr(offdesign=["m"])
#gt1.conns["c3"].set_attr(design=["T"]) # because I want variable TIT
c18.set_attr(design=["p"])
hs.set_attr(design=["pr2"], offdesign=["zeta1", "zeta2"]) # "kA_char" ?
pipe_feed.set_attr(design=["pr"], offdesign=["zeta"]) # "Q", "kA_char", kA_group ?
cons.set_attr(design=["pr"], offdesign=["zeta"])
#val.set_attr(offdesign=["zeta"])
pipe_return.set_attr(design=["pr"], offdesign=["zeta"]) # "Q", "kA_char", kA_group ?

print("Off-design performance")
nw.solve("offdesign", design_path="CHP & DHN_identical GTs_design point")
nw.print_results()
print(f"Brayton cycle 1 net electric power output = {round(abs(generator1.P.val)/1e6, 3)} MWe", "\n")
el_eff1 = (abs(generator1.P.val) / gt1.comps["combustion chamber"].ti.val) * 100 # Brayton cycle (electric) efficiency
print(f"Brayton cycle 1 (electric) efficiency = {round(el_eff1, 2)} %", "\n")
print(c10.fluid.val)
print(c17.fluid.val)
print(f"DHN mass flow rate = {round(c0.m.val, 3)} kg/s")
print(f"Thermal power given from CHP plant to (primary) DHN = {round(abs(hs.Q.val)/1e6, 3)} MWth")
print(f"Thermal power given from primary DHN to secondary DHN (consumers) = {round(abs(cons.Q.val)/1e6, 3)} MWth")
print(f"Thermal efficiency of (primary) DHN = {round((abs(cons.Q.val)/abs(hs.Q.val)) * 100, 2)} %")
print("\n")

nw.set_attr(iterinfo=False)

p_ambient = 0.93194 # [bar]
T_ambient_list = [13.5, 8, 3.9, 2.3, 3.7, 6.9, 11.6, 16.8] # [degrees C]
# Gas turbines' combined net electric power output
GTs_power = []
# Thermal power given from CHP plant to (primary) DHN
CHP_to_DHN = []
n = 1

for T in T_ambient_list:
    T_ambient = T
    gt1.conns["c1"].set_attr(p=p_ambient, T=T_ambient)
    print(f"Month: {n}")
    nw.solve("offdesign", design_path="CHP & DHN_identical GTs_design point")
    nw.print_results()
    GTs_power += [(2*abs(generator1.P.val)) / 1e6] # [MWe]
    CHP_to_DHN += [abs(hs.Q.val) / 1e6] # [MWth]
    print(f"Brayton cycle 1 net electric power output = {round(abs(generator1.P.val)/1e6, 3)} MWe", "\n")
    el_eff1 = (abs(generator1.P.val) / gt1.comps["combustion chamber"].ti.val) * 100
                                                # Brayton cycle (electric) efficiency
    print(f"Brayton cycle 1 (electric) efficiency = {round(el_eff1, 2)} %", "\n")
    print(c10.fluid.val)
    print(c17.fluid.val)
    print(f"DHN mass flow rate = {round(c0.m.val, 3)} kg/s")
    print(f"Thermal power given from CHP plant to (primary) DHN = {round(abs(hs.Q.val) / 1e6, 3)} MWth")
    print(f"Thermal power given from primary DHN to secondary DHN (consumers) = {round(abs(cons.Q.val) / 1e6, 3)} MWth")
    print(f"Thermal efficiency of (primary) DHN = {round((abs(cons.Q.val) / abs(hs.Q.val)) * 100, 2)} %")
    print("\n")
    n += 1

# plotting
import matplotlib.pyplot as plt

fig, ax = plt.subplots(2, 1, figsize=(16, 8), sharex='col')

ax = ax.flatten()
[a.grid() for a in ax]

ax[0].scatter(T_ambient_list, GTs_power, s=100, color="#1f567d")
ax[1].scatter(T_ambient_list, CHP_to_DHN, s=100, color="#18a999")

ax[0].set_ylabel('CHP plant net electric power output [MWe]')
ax[1].set_ylabel('Thermal power given from CHP plant to DHN [MWth]')
ax[1].set_xlabel('Ambient temperature [oC]')

plt.tight_layout()
fig.savefig('Ambient temperature_off-design.svg')
plt.close()

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

CHP plant & DHN design_Problem setting up the system #476

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 2 comments 5 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

Select a reply

CHP plant & DHN design_Problem setting up the system #476

MikeSennis Jan 17, 2024

Replies: 2 comments · 5 replies

fwitte Jan 21, 2024 Maintainer

fwitte Jan 21, 2024 Maintainer

fwitte Jan 21, 2024 Maintainer

MikeSennis Jan 24, 2024 Author

fwitte Jan 24, 2024 Maintainer

MikeSennis Jan 25, 2024 Author

MikeSennis Feb 5, 2024 Author

MikeSennis
Jan 17, 2024

Replies: 2 comments 5 replies

fwitte
Jan 21, 2024
Maintainer

fwitte Jan 21, 2024
Maintainer

fwitte Jan 21, 2024
Maintainer

MikeSennis
Jan 24, 2024
Author

fwitte Jan 24, 2024
Maintainer

MikeSennis Jan 25, 2024
Author

MikeSennis Feb 5, 2024
Author