oemof / tespy

Thermal Engineering Systems in Python (TESPy). This package provides a powerful simulation toolkit for thermal engineering plants such as power plants, district heating systems or heat pumps.
https://tespy.readthedocs.io
MIT License
272 stars 85 forks source link

Drum component fluid value error: parallel heat pumps #274

Closed NicholasFry closed 3 years ago

NicholasFry commented 3 years ago

Hi there. I am modifying the COP of a heat pump tutorial to create parallel heat pumps feeding a district heating network, using a waste heat source of the same rejection reservoir. It seems if I remove air from the tutorial, the network solver works fine. If I remove a compressor, an intercooler, and air as a fluid from the code, then I get a problem with the evaporator drum.

ValueError: The function h_mix_pQ can only be used for pure fluids.

I am not certain that I understand the conflict in this case. My fractions of NH3 and water are not mixed at any inlets. Hopefully this is the appropriate place to raise the question.


#the goal for this will be two heat pumps in "supercharged configuration" where stage 2 starts the train from the cooling reservoir
# -*- coding: utf-8 -*-
from tespy.networks import Network
from tespy.components import (
    Sink, Source, Splitter, Compressor, Condenser, Pump, HeatExchangerSimple,
    Valve, Drum, HeatExchanger, CycleCloser
)
from tespy.connections import Connection, Ref
from tespy.tools.characteristics import CharLine
from tespy.tools.characteristics import load_default_char as ldc
from tespy.tools import document_model

import numpy as np
import pandas as pd

# %% network
#Water is used for the cold side of the heat exchanger, for the consumer and for the hot side of the environmental temperature.
#Ammonia is used as coolant within the heat pump circuit.
nw = Network(
    fluids=['water', 'NH3'], T_unit='C', p_unit='bar', h_unit='kJ / kg',
    m_unit='kg / s'
)

# %% components

# sources & sinks
cc1 = CycleCloser('coolant cycle closer1')#CycleCloser component makes sure, the fluid properties pressure and enthalpy are identical at the inlet and the outlet. 
cc2 = CycleCloser('coolant cycle closer2')
cc_cons = CycleCloser('consumer cycle closer')
electra_therm_rejected_heat = Source('source heat from electratherm')
# exit_heat_from_condenser1 = Source('source heat from condenser1 exit')#added by trial and error
# exit_heat_from_evaporator1 = Sink('exit_heat_from_evaporator1')#added by trial and error
sink_electratherm_cooling_reservoir= Sink('sink cool water')

# ambient system
sp = Splitter('splitter')#I should not need this because it was originally to split between the intercooler and the evaporator supply
pu = Pump('pump')#pumps source heat from ElectraTherm rejected heat reservoir

# consumer system

condenser_1 = Condenser('condenser 1')#closest to cooling reservoir
condenser_2 = Condenser('condenser 2')
district_heating_pump = Pump('district heating pump')#feeds the first condenser 
cons_1 = HeatExchangerSimple('consumer 1')#exits the first condenser, enters the second condenser
# cons_2 = HeatExchangerSimple('consumer 2')#exits the second condenser, supply temperature for DH network
#I may not need cons_2 because of the consumer fluid connection sequence
# evaporator system(s)

ves1 = Valve('valve 1')
ves2 = Valve('valve 2')
dr1 = Drum('drum 1')
dr2 = Drum('drum 2')
ev1 = HeatExchanger('evaporator 1')
ev2 = HeatExchanger('evaporator 2')
su1 = HeatExchanger('superheater 1')
su2 = HeatExchanger('superheater 2')
erp1 = Pump('evaporator reciculation pump 1')
erp2 = Pump('evaporator reciculation pump 2')

# compressor-system(s)

compressor1 = Compressor('compressor 1')
compressor2 = Compressor('compressor 2')
# compressor3 = Compressor('compressor 3')
# compressor4 = Compressor('compressor 4')
# ic1 = HeatExchanger('intercooler 1')#I have removed the intercooler, an intermediate step between two compressors in a single unit
# ic2 = HeatExchanger('intercooler 2')

# %% connections

# consumer system

c_in_condenser_1 = Connection(cc1, 'out1', condenser_1, 'in1')
c_in_condenser_2 = Connection(cc2, 'out1', condenser_2, 'in1')

cb_district_heating_pump = Connection(cc_cons, 'out1', district_heating_pump, 'in1')#return from DH
district_heating_pump_condenser_1 = Connection(district_heating_pump, 'out1', condenser_1, 'in2')#DH pump feeds condenser 1 (return from DH)
cons_condenser1_feeds_condenser2 = Connection(condenser_1, 'out2', condenser_2, 'in2')
condenser2_cons = Connection(condenser_2, 'out2', cons_1, 'in1')#condenser_2 out to DH
# condenser2_cons = Connection(condenser_2, 'out2', cons_2, 'in1')#condenser_2 out to consumer
# cons1_out_to_cons2_in = Connection(cons_1, 'out1', cons_2, 'in1')#consumer side between condensers
cons_cf = Connection(cons_1, 'out1', cc_cons, 'in1')#out to consumer DH network with cycle closer for enthalpy and pressure match

nw.add_conns(c_in_condenser_1, c_in_condenser_2, cb_district_heating_pump, district_heating_pump_condenser_1, cons_condenser1_feeds_condenser2, condenser2_cons, cons_cf)

# connection condenser - evaporator system

condenser1_ves = Connection(condenser_1, 'out1', ves1, 'in1')#condeser to expansion valve
condenser2_ves = Connection(condenser_2, 'out1', ves2, 'in1')

nw.add_conns(condenser1_ves, condenser2_ves)

# evaporator system

ves1_dr1 = Connection(ves1, 'out1', dr1, 'in1')#valve connects to drum at inlet
ves2_dr2 = Connection(ves2, 'out1', dr2, 'in1')#valve connects to drum at inlet
dr1_erp1 = Connection(dr1, 'out1', erp1, 'in1')#drum outlet to evaporator recirculation pump inlet
dr2_erp2 = Connection(dr2, 'out1', erp2, 'in1')#drum outlet to evaporator recirculation pump inlet
erp1_ev1 = Connection(erp1, 'out1', ev1, 'in2')#pump outlet to evaporator inlet
erp2_ev2 = Connection(erp2, 'out1', ev2, 'in2')#pump outlet to evaporator inlet
ev1_dr1 = Connection(ev1, 'out2', dr1, 'in2')#evaporator outlet to drum inlet
ev2_dr2 = Connection(ev2, 'out2', dr2, 'in2')#evaporator outlet to drum inlet
dr1_su1 = Connection(dr1, 'out2', su1, 'in2')#drum outlet to the 
dr2_su2 = Connection(dr2, 'out2', su2, 'in2')#drum outlet to the superheater2 (su2 is not currently doing anything because we will not connect another pump to it, I think)

nw.add_conns(ves1_dr1, ves2_dr2, dr1_erp1, dr2_erp2, erp1_ev1, erp2_ev2, ev1_dr1, ev2_dr2, dr1_su1, dr2_su2)#add all the connections above to the network, Connections is a class, add_conns is the function that takes them in

rejected_heat_to_pump = Connection(electra_therm_rejected_heat, 'out1', pu, 'in1')
p_su1 = Connection(pu, 'out1', su1, 'in1')#pump to superheater1
# sp_su = Connection(sp, 'out1', su, 'in1')
su1_ev1 = Connection(su1, 'out1', ev1, 'in1')#superheater1 to evaporator1
ev1_su2 = Connection(ev1, 'out1', su2, 'in1')
su2_ev2 = Connection(su2, 'out1', ev2, 'in1')
# ev1_ev2 = Connection(ev1, 'out1', ev2, 'in1')#connect ev1 heat rejection to ev2 inlet
# ev1_sink = Connection(ev1, 'out1', sink_electratherm_cooling_reservoir, 'in1')#not sure on this yet
ev2_sink = Connection(ev2, 'out1', sink_electratherm_cooling_reservoir, 'in1')#Does this work????

nw.add_conns(rejected_heat_to_pump, p_su1, su1_ev1, su2_ev2, ev1_su2, ev2_sink)

# connection evaporator system - compressor system

su1_compressor1 = Connection(su1, 'out2', compressor1, 'in1')
su2_compressor2 = Connection(su2, 'out2', compressor2, 'in1')#recall that su2 is not connected to any heat source at the moment

nw.add_conns(su1_compressor1, su2_compressor2)

# compressor-system

# compressor1_he = Connection(compressor1, 'out1', ic, 'in1')#compressor1 heat exchanger, comment out if no intercooler or twin compressor train
# he_compressor2 = Connection(ic, 'out1', compressor2, 'in1')#compressor2 heat exchanger, comment out if no intercooler or twin compressor train
compressor1_c_out = Connection(compressor1, 'out1', cc1, 'in1')
compressor2_c_out = Connection(compressor2, 'out1', cc2, 'in1')

# sp_ic = Connection(sp, 'out2', ic, 'in2')
# ic_out = Connection(ic, 'out2', amb_out2, 'in1')

nw.add_conns(compressor1_c_out, compressor2_c_out)

# %% component parametrization

# condenser system

condenser_1.set_attr(pr1=0.99, pr2=0.99, ttd_u=5, design=['pr2', 'ttd_u'], #upper terminal temperature difference as design parameter, pressure ratios
            offdesign=['zeta2', 'kA_char'])#kA_char is area independent heat transfer coefficient characteristic
condenser_2.set_attr(pr1=0.99, pr2=0.99, ttd_u=5, design=['pr2', 'ttd_u'], #upper terminal temperature difference as design parameter, pressure ratios
            offdesign=['zeta2', 'kA_char'])#kA_char is area independent heat transfer coefficient characteristic
district_heating_pump.set_attr(eta_s=0.8, design=['eta_s'], offdesign=['eta_s_char'])#efficiency of the district heating pump set to 80%
cons_1.set_attr(pr=0.99, design=['pr'], offdesign=['zeta'])#In offdesign calculation the consumer’s pressure ratio will be a function of the mass flow, thus as offdesign parameter we select zeta
# cons_2.set_attr(pr=0.99, design=['pr'], offdesign=['zeta'])#In offdesign calculation the consumer’s pressure ratio will be a function of the mass flow, thus as offdesign parameter we select zeta

# water pump

pu.set_attr(eta_s=0.75, design=['eta_s'], offdesign=['eta_s_char'])#this is for a 75% pump efficiency

# evaporator system

kA_char1 = ldc('heat exchanger', 'kA_char1', 'DEFAULT', CharLine)#Characteristic line for hot side heat transfer coefficient. heat transfer coefficient multiplied by the area of HEX, Charline is the linear interpolation equation of the x, y; DEFAULT HEX function is charted here https://tespy.readthedocs.io/en/main/api/tespy.data.html?highlight=DEFAULT#default-characteristics
kA_char2 = ldc('heat exchanger', 'kA_char2', 'EVAPORATING FLUID', CharLine)#Characteristic line for cold side heat transfer coefficient (ammonia in this case for Mandaree)

ev1.set_attr(pr1=0.98, pr2=0.99, ttd_l=5, #pressure_ratio1, pressure_ratio2, terminal_temperature_difference1(terminal temperature difference at the evaporator’s cold side inlet)
            kA_char1=kA_char1, kA_char2=kA_char2,
            design=['pr1', 'ttd_l'], offdesign=['zeta1', 'kA_char'])
ev2.set_attr(pr1=0.98, pr2=0.99, ttd_l=5, #pressure_ratio1, pressure_ratio2, terminal_temperature_difference1(terminal temperature difference at the evaporator’s cold side inlet)
            kA_char1=kA_char1, kA_char2=kA_char2,
            design=['pr1', 'ttd_l'], offdesign=['zeta1', 'kA_char'])
su1.set_attr(pr1=0.98, pr2=0.99, ttd_u=2, design=['pr1', 'pr2', 'ttd_u'], #ttd_u refers to the boiling temperature at hot side inlet
            offdesign=['zeta1', 'zeta2', 'kA_char'])
su2.set_attr(pr1=0.98, pr2=0.99, ttd_u=2, design=['pr1', 'pr2', 'ttd_u'],
            offdesign=['zeta1', 'zeta2', 'kA_char'])
erp1.set_attr(eta_s=0.8, design=['eta_s'], offdesign=['eta_s_char'])#eta_s is the Isentropic (adiabatic) efficiency of the process
erp2.set_attr(eta_s=0.8, design=['eta_s'], offdesign=['eta_s_char'])

# compressor system

compressor1.set_attr(eta_s=0.85, design=['eta_s'], offdesign=['eta_s_char'])
compressor2.set_attr(eta_s=0.85, pr=3, design=['eta_s'], offdesign=['eta_s_char'])#pressure ratio 3:1 outlet:inlet
# ic.set_attr(pr1=0.99, pr2=0.98, design=['pr1', 'pr2'],
#             offdesign=['zeta1', 'zeta2', 'kA_char'])

# %% connection parametrization

# condenser system

c_in_condenser_1.set_attr(fluid={'NH3': 1, 'water': 0})
cb_district_heating_pump.set_attr(T=60, p=10, fluid={'NH3': 0, 'water': 1})
condenser2_cons.set_attr(T=90)

# evaporator system cold side

erp1_ev1.set_attr(m=Ref(ves1_dr1, 1.25, 0), p0=5)#pump outlet to evaporator1
erp2_ev2.set_attr(m=Ref(ves2_dr2, 1.25, 0), p0=5)#pump outlet to evaporator2
su1_compressor1.set_attr(p0=5, state='g')
su2_compressor2.set_attr(p0=5, state='g')#p0 starting value specification for pressure.

# evaporator system hot side (from ElectraTherm rejection)

# pumping at constant rate in partial load
rejected_heat_to_pump.set_attr(T=12, p=2, fluid={'NH3': 0, 'water': 1},
               offdesign=['v'])
# sp_su.set_attr(offdesign=['v'])
ev2_sink.set_attr(p=2, T=9, design=['T'])

# # compressor-system intercooler inlet and outlet; commenting out 

# he_compressor2.set_attr(Td_bp=5, p0=20, design=['Td_bp'])
# ic_out.set_attr(T=30, design=['T'])

# %% key paramter (consumer demand)

cons_1.set_attr(Q=6000000) #6MW demand, demand unit in watts

# %% Calculation

nw.solve('design')#network solve    
nw.print_results()
nw.save('heat_pump_water')
document_model(nw, filename='report_water_design.tex')#output network model to latex report

# offdesign test
nw.solve('offdesign', design_path='heat_pump_water')#solve the offdesign values for the network (other projected outcomes)
document_model(nw, filename='report_water_offdesign.tex')#print these alternatives to a latex report

T_range = [6, 12, 18, 24, 30]
Q_range = np.array([100e3, 120e3, 140e3, 160e3, 180e3, 200e3, 220e3])
df = pd.DataFrame(columns=Q_range / -cons_1.Q.val)

for T in T_range:
    rejected_heat_to_pump.set_attr(T=T)
    eps = []

    for Q in Q_range:
        cons_1.set_attr(Q=-Q)
        nw.solve('offdesign', design_path='heat_pump_water')

        if nw.lin_dep:
            eps += [np.nan]
        else:
            eps += [
                abs(condenser_1.Q.val+condenser_2.Q.val) / (compressor1.P.val + compressor2.P.val + erp1.P.val + erp2.P.val + pu.P.val)
            ]

    df.loc[T] = eps

df.to_csv('COP_water.csv')
fwitte commented 3 years ago

Hi @NicholasFry,

you have missed to specify the fluid in the second heat pump cycle. That seems to be the reason for the error message. However, I think there should be a different error message in this case. Basically, the h_mix_pQ function is called to generate starting values for the solver while missing the fluid information. If you add

c_in_condenser_2.set_attr(fluid={'NH3': 1, 'water': 0})

at least that issue is fixed, but the solver does not find a solution. Might be due to bad specifications or due to bad starting values. As the network is quite large it is always difficult to search the error if you did not set up the network yourself, so I cannot immediately tell you what the issue is. I'd advise you to set up your model step by step. Use one heat pump only and the add component by component to the secondary cycle, i.e. insert the second condenser inside the district heating but only connect it to a source and a sink on the heat pump side. Then add the valve, the evaporator, ... :).

If you need any more help or support, just let me know :)

Best regards, have a nice week Francesco

P.S.: That looks like an interesting set up! Analyzing the influence of using different working fluids for the two heat pumps on the overall performance could be very interesting.

NicholasFry commented 3 years ago

Thank you for identifying the next steps for me @fwitte. I am enjoying learning about the capabilities of TESPy. I appreciate the feedback.