oemof / oemof-examples

A collection of oemof examples and notebooks.
https://oemof.org/
MIT License
44 stars 69 forks source link

Starting parameter assistance #76

Open NicholasFry opened 3 years ago

NicholasFry commented 3 years ago

Hi all. I am looking for more assistance on defining the initial parameters in this basic heat pump layout. This is fed by hot water from a waste heat source.

What do I need to do to define the appropriate starting values (e.g. pr=, ttd_u=, p0= outlet to the evaporator)? Can I incorporate the outdoor air temperature conditions?

How can I determine the correct mass flow factor below, prior to running the calculations? erp1_ev1.set_attr(m=Ref(ves1_dr1, 1.25, 0), p0=5)#pump outlet to evaporator1

With a 4MWthermal design size, how can I find an optimum if the COP keeps increasing with smaller relative load factors?

Thank you for any assistance you can offer. I am still trying to learn how to appropriately design with this library.

# -*- 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')#add later for heat pump 2
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')
pu = Pump('pump')#pumps source heat from ElectraTherm rejected heat reservoir

# consumer system

condenser_1 = Condenser('condenser 1')#closest to cooling reservoir
district_heating_pump = Pump('district heating pump')#feeds the first condenser 
cons_1 = HeatExchangerSimple('consumer 1')#exits the first condenser, enters the second condenser

ves1 = Valve('valve 1')
dr1 = Drum('drum 1')
ev1 = HeatExchanger('evaporator 1')
erp1 = Pump('evaporator recirculation pump 1')

# compressor-system(s)

compressor1 = Compressor('compressor 1')

# %% connections

# consumer system

c_in_condenser_1 = Connection(cc1, 'out1', condenser_1, 'in1')#cyclecloser to condenser 1, see https://tespy.readthedocs.io/en/main/tutorials_examples.html#cop-of-a-heat-pump

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)
condenser1_cons = Connection(condenser_1, 'out2', cons_1, 'in1')
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, cb_district_heating_pump, district_heating_pump_condenser_1, condenser1_cons, cons_cf)

# connection condenser - evaporator system

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

nw.add_conns(condenser1_ves)

ves1_dr1 = Connection(ves1, 'out1', dr1, 'in1')#valve connects to drum at inlet
dr1_erp1 = Connection(dr1, 'out1', erp1, 'in1')#drum outlet to evaporator recirculation pump inlet
erp1_ev1 = Connection(erp1, 'out1', ev1, 'in2')#pump outlet to evaporator inlet
ev1_dr1 = Connection(ev1, 'out2', dr1, 'in2')#evaporator outlet to drum inlet

nw.add_conns(ves1_dr1, dr1_erp1, erp1_ev1, ev1_dr1)#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_ev1 = Connection(pu, 'out1', ev1, 'in1')#pump to evaporator1
ev1_sink = Connection(ev1, 'out1', sink_electratherm_cooling_reservoir, 'in1')#exit to cooling reservoir

nw.add_conns(rejected_heat_to_pump, p_ev1, ev1_sink)

# connection evaporator system - compressor system

dr1_compressor1 = Connection(dr1, 'out2', compressor1, 'in1')#drum outlet2 to compressor1

nw.add_conns(dr1_compressor1)

compressor1_c_out = Connection(compressor1, 'out1', cc1, 'in1')

nw.add_conns(compressor1_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
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

# 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'])
erp1.set_attr(eta_s=0.8, design=['eta_s'], offdesign=['eta_s_char'])#eta_s is the Isentropic (adiabatic) efficiency of the process

# compressor system

compressor1.set_attr(eta_s=0.85, design=['eta_s'], offdesign=['eta_s_char'])#docs say not to set compressor 1 pressure ratio for parallel, if adding another compressor later remove pressure ratio, #pressure ratio 3:1 outlet:inlet

# %% connection parametrization

# condenser system

c_in_condenser_1.set_attr(fluid={'NH3': 1, 'water': 0})
cb_district_heating_pump.set_attr(T=50, p=10, fluid={'NH3': 0, 'water': 1})
condenser1_cons.set_attr(T=70)#tempearture feeding the DH network

# evaporator system cold side

erp1_ev1.set_attr(m=Ref(ves1_dr1, 1.25, 0), p0=5)#pump outlet to evaporator1 Ref(ref_obj, factor, delta)

# evaporator system hot side (from ElectraTherm rejection)

# pumping at constant rate in partial load
rejected_heat_to_pump.set_attr(T=60, p=2, fluid={'NH3': 0, 'water': 1},
               offdesign=['v'])#here I assume we are maintaining 2 bar in the partial load condition from rejected heat reservoir from generators at 60C
ev1_sink.set_attr(p=2, T=20, design=['T'])

# %% key paramter (consumer demand)

cons_1.set_attr(Q=-4000000) #4MW demand, demand unit in watts, this value should be negative to represent a demand on the system

# %% 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 = [35, 40, 45, 50, 55, 65]
Q_range = np.array([50e4, 180e4, 280e4, 300e4, 340e4, 360e4, 380e4, 420e4])
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) / (compressor1.P.val + erp1.P.val + pu.P.val)
            ]

    df.loc[T] = eps

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

Hi @NicholasFry, I will be on a short train ride later and have a look at the issue. Best Francesco

fwitte commented 3 years ago

Btw: I saw that you are interested in geothermal heating. We have been working on geothermal power generation recently (https://github.com/fwitte/ORCSimulator) as well as geothermal storage integration with heat pumps. If you are interested, we can have a talk on that topic. Just contact me via matrix (https://matrix.to/#/@fwitte:matrix.org) or email (francesco.witte( at )hs-flensburg.de). Maybe there might be some idea for a cooperation?

fwitte commented 3 years ago

What do I need to do to define the appropriate starting values (e.g. pr=, ttd_u=, p0= outlet to the evaporator)?

Starting values are only required for the Newton-Algorithm as good initial guess for faster (often also only way of successful) solution finding. The model you apply does not have starting value issues in setting up the design but rather than that in the offdesign simulation. As we have a good initial guess from the design simulation, states, that are too far away from that point might in some cases be difficult to find a solution for. I recommend a procedure similar to the following:

[...]
T_range = [35, 40, 45, 50, 55, 65][::-1]
Q_range = np.array([50e4, 180e4, 280e4, 300e4, 340e4, 360e4, 380e4, 420e4])[::-1]
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)
        if Q == Q_range[0]:
            nw.solve('offdesign', design_path='heat_pump_water', init_path='heat_pump_water')
        else:
            nw.solve('offdesign', design_path='heat_pump_water')

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

    df.loc[T] = eps

df.to_csv('COP_water.csv')

Can I incorporate the outdoor air temperature conditions?

What do you mean here? For what purpose do you need the ambient air conditions

How can I determine the correct mass flow factor below, prior to running the calculations? erp1_ev1.set_attr(m=Ref(ves1_dr1, 1.25, 0), p0=5)#pump outlet to evaporator1

This value basically does nothing in this model besides determining the steam mass fraction of the working fluid returning from the evaporator into the drum. I kind of set this value arbitrarily. I actually not even sure, if using a drum for gas/liquid separation is commong in heat pumps.

With a 4MWthermal design size, how can I find an optimum if the COP keeps increasing with smaller relative load factors?

TESPy does not have feature inbuilt for this function. But you can easily include the simulation run in some kind of search method like the golden section search or bisection. Have a look into the ORC github repository for an implementation of the golden section search. The code is deriven from the example code at https://en.wikipedia.org/wiki/Golden-section_search.

NicholasFry commented 3 years ago

@fwitte thank you again for answering my questions. I will look into the optimization process. As for the drum separator, I do not have the answer on typical designs either. I appreciate the insights on how it works.

My ambient temperature concern is more of a feeding fluid temperature problem. Since I am drawing waste heat from a fluid reservoir, it will be subject to the air temperature. In this case, I should simply vary the feed temperature profile with climatic conditions.

As for the 35C case, it then seems possible to use TESPy to quickly determine the lower limit for feed temperature on a cycle, given the stated exit temperature. This is very interesting to work with and actually quite accessible! Enjoy your weekend.