no simple heat pump model solution when changing fluid from NH3 to CO2 #352

matlaw · 2022-08-14T05:19:48Z

matlaw
Aug 14, 2022

Using the fluprodia example file for plotting states into a fluid property diagram found here: https://fluprodia.readthedocs.io/en/latest/usage.htm, everything works for fluid NH3 (see working code below).

However, when I change the two references to NH3 to CO2, I get this error:

ValueError: unable to solve 1phase PY flash with Tmin=216.592, Tmax=3000 due to error: HSU_P_flash_singlephase_Brent could not find a solution because Hmolar [13693.8 J/mol] is below the minimum value of 19160.774415 J/mol

I'm assuming it has to do with lines:
cd_va.set_attr(Td_bp=-5, T=85)
ev_cp.set_attr(Td_bp=5, T=15)

but I'm not really sure what values I should be changing them to when switching fluid.

Any suggestions would be appreciated - thanks!

from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Valve)
from tespy.connections import Connection
from tespy.networks import Network

def run_simple_heat_pump_model():
    nw = Network(['NH3'], T_unit='C', p_unit='bar', h_unit='kJ / kg')
    nw.set_attr(iterinfo=False)
    cp = Compressor('compressor')
    cc = CycleCloser('cycle_closer')
    cd = HeatExchangerSimple('condenser')
    va = Valve('expansion valve')
    ev = HeatExchangerSimple('evaporator')

    cc_cd = Connection(cc, 'out1', cd, 'in1')
    cd_va = Connection(cd, 'out1', va, 'in1')
    va_ev = Connection(va, 'out1', ev, 'in1')
    ev_cp = Connection(ev, 'out1', cp, 'in1')
    cp_cc = Connection(cp, 'out1', cc, 'in1')

    nw.add_conns(cc_cd, cd_va, va_ev, ev_cp, cp_cc)

    cd.set_attr(pr=0.95, Q=-1e6)
    ev.set_attr(pr=0.9)
    cp.set_attr(eta_s=0.9)

    cc_cd.set_attr(fluid={'NH3': 1})
    cd_va.set_attr(Td_bp=-5, T=85)
    ev_cp.set_attr(Td_bp=5, T=15)

    nw.solve('design')

    result_dict = {}
    result_dict.update(
        {cp.label: cp.get_plotting_data()[1] for cp in nw.comps['object']
         if cp.get_plotting_data() is not None})

    return result_dict

tespy_results = run_simple_heat_pump_model()

fwitte · 2022-08-14T06:49:21Z

fwitte
Aug 14, 2022
Maintainer

Good morning @matlaw,

thank you for reaching out. You are looking at the correct place, the issue is, that the CO2 is mostly at supercritical tempertures given the inputs (see http://coolprop.org/fluid_properties/fluids/CarbonDioxide.html#fluid-carbondioxide). So you cannot specify the temperature difference to the boiling point, because it does not exist after "condensation". If you want to model a heat pump with (partially) supercrticial working fluid, you'd need to make some small adjustments to the setup. I attached an example below, that does provide a stable start for either a supercritical or a subcritical cycle. I am not yet very familiar with supercritical heat pumps, so I do not really know whether my specifications would make sense like this in an actually desigining the plant.

In generaly, note

Remove the isentropic efficiency specification in the setup of the plant, it often causes trouble
Apply it back after getting the stable starting values
Use pressure and enthalpy specifications for generating good starting values, because these are the primariy variables of the solver (that means, the equations are extremly easily solved)
The temperature levels applied do not really make sense. The power input of the compressor is negative because the cycle is actually inverted. The enthalpy after the heat discharge is higher than the enthalpy after the compression.

from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Valve)
from tespy.connections import Connection
from tespy.networks import Network
from CoolProp.CoolProp import PropsSI as PSI


def run_simple_heat_pump_model(working_fluid, supercritical, T_hot, T_cold):
    nw = Network([working_fluid], T_unit='C', p_unit='bar', h_unit='kJ / kg')
    nw.set_attr(iterinfo=False)
    cp = Compressor('compressor')
    cc = CycleCloser('cycle_closer')
    cd = HeatExchangerSimple('condenser')
    va = Valve('expansion valve')
    ev = HeatExchangerSimple('evaporator')

    cc_cd = Connection(cc, 'out1', cd, 'in1', label='0')
    cd_va = Connection(cd, 'out1', va, 'in1', label='1')
    va_ev = Connection(va, 'out1', ev, 'in1', label='2')
    ev_cp = Connection(ev, 'out1', cp, 'in1', label='3')
    cp_cc = Connection(cp, 'out1', cc, 'in1', label='4')

    nw.add_conns(cc_cd, cd_va, va_ev, ev_cp, cp_cc)

    cd.set_attr(pr=1, Q=-1e6)
    ev.set_attr(pr=1)

    cc_cd.set_attr(fluid={working_fluid: 1})

    if supercritical:
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot > T_critical, "Hot side temperature specification is below critical temperature"
        assert T_cold > T_critical, "Cold side temperature specification is below critical temperature"

        # for a good guess of pressure levels
        p_critical = PSI('p_critical', working_fluid)
        p_high = p_critical * 2
        # you do not want to be exactly at the critical point, here CoolProp ofte has difficulties solving
        # the equation of state
        p_low = p_critical * 1.1

        # for a good guess of enthalpy values
        h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid)
        h_after_heat_input = PSI('H', 'P', p_low, 'T', T_cold + 273.15, working_fluid)
        h_after_compression = h_after_heat_input + 50e3

        cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3)
        ev_cp.set_attr(p=p_low / 1e5, h=h_after_heat_input / 1e3)
        cp_cc.set_attr(h=h_after_compression / 1e3)
    else:
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot < T_critical, "Hot side temperature specification is above critical temperature"
        assert T_cold < T_critical, "Cold side temperature specification is above critical temperature"

        h_condensate = PSI('H', 'Q', 0, 'T', T_hot + 273.15, working_fluid)
        p_condensate = PSI('P', 'Q', 0, 'T', T_hot + 273.15, working_fluid)
        h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        h_after_compression = h_evaporation + 50e3
        cd_va.set_attr(p=p_condensate / 1e5, h=h_condensate / 1e3)
        ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3)
        cp_cc.set_attr(h=h_after_compression / 1e3)

    nw.solve('design')

    return nw


def extract_plotting_results(nw):

    result_dict = {}
    result_dict.update(
        {cp.label: cp.get_plotting_data()[1] for cp in nw.comps['object']
         if cp.get_plotting_data() is not None})

    return result_dict


tespy_stable_nw = run_simple_heat_pump_model('CO2', True, 85, 40)
c1, c3, c4 = tespy_stable_nw.get_conn(['1', '3', '4'])
comp = tespy_stable_nw.get_comp('compressor')

# unset stable start value specifications and set target values instead
# for example compressor efficiency
c4.set_attr(h=None)
comp.set_attr(eta_s=0.9)

# or actual temperature levels instead of enthalpy
c1.set_attr(h=None, T=85)
c3.set_attr(h=None, T=40)
tespy_stable_nw.solve('design')

tespy_stable_nw.print_results()

Generally, the solving algorithm is quite senstive to a good initial guess in TESPy. A more extensive example on setting up a subcritical heat pump with easily interchangable working fluids had been started at #346, I'll make an update of the TESPy online documentation soon to include that.

Maybe we should add supercritical and transcritical to that as well. If you plan to work more with TESPy in the future and are interested in contributing you are welcome to reach out anytime :).

Have a nice sunday!

3 replies

matlaw Aug 15, 2022
Author

Thanks a lot Francesco, that was really helpful.

Could you point me in the direction of what I would do to enable modelling of a transcritical cycle?

Matthew

fwitte Aug 15, 2022
Maintainer

Great! I‘ll write the something up for you later today:)

fwitte Aug 17, 2022
Maintainer

Sorry for the lengthy wait @matlaw, here you go:

For transcritical case you want to use the same specifications as supercritical for heat discharge and as subcritical for evaporation
When changing up parameters given the stable starting values, set fully saturated gas after evaporation and desired temperature level of evaporation instead of pressure and enthalpy

from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Valve)
from tespy.connections import Connection
from tespy.networks import Network
from CoolProp.CoolProp import PropsSI as PSI


def run_simple_heat_pump_model(working_fluid, mode, T_hot, T_cold):
    nw = Network([working_fluid], T_unit='C', p_unit='bar', h_unit='kJ / kg')
    nw.set_attr(iterinfo=False)
    cp = Compressor('compressor')
    cc = CycleCloser('cycle_closer')
    cd = HeatExchangerSimple('condenser')
    va = Valve('expansion valve')
    ev = HeatExchangerSimple('evaporator')

    cc_cd = Connection(cc, 'out1', cd, 'in1', label='0')
    cd_va = Connection(cd, 'out1', va, 'in1', label='1')
    va_ev = Connection(va, 'out1', ev, 'in1', label='2')
    ev_cp = Connection(ev, 'out1', cp, 'in1', label='3')
    cp_cc = Connection(cp, 'out1', cc, 'in1', label='4')

    nw.add_conns(cc_cd, cd_va, va_ev, ev_cp, cp_cc)

    cd.set_attr(pr=1, Q=-1e6)
    ev.set_attr(pr=1)

    cc_cd.set_attr(fluid={working_fluid: 1})

    if mode == "supercritical":
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot > T_critical, "Hot side temperature specification needs to be above critical temperature"
        assert T_cold > T_critical, "Cold side temperature specification needs to be above critical temperature"

        # for a good guess of pressure levels
        p_critical = PSI('p_critical', working_fluid)
        p_high = p_critical * 2
        # you do not want to be exactly at the critical point, here CoolProp ofte has difficulties solving
        # the equation of state
        p_low = p_critical * 1.1

        # for a good guess of enthalpy values
        h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid)
        h_after_heat_input = PSI('H', 'P', p_low, 'T', T_cold + 273.15, working_fluid)

        cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3)
        ev_cp.set_attr(p=p_low / 1e5, h=h_after_heat_input / 1e3)
    elif mode == "transcritical":
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot > T_critical, "Hot side temperature specification needs to be above critical temperature"
        assert T_cold < T_critical, "Cold side temperature specification needs to be below critical temperature"

        # for a good guess of pressure levels
        p_critical = PSI('p_critical', working_fluid)
        p_high = p_critical * 2
        h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid)
        h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3)
        ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3)

    else:
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot < T_critical, "Hot side temperature specification needs to be below critical temperature"
        assert T_cold < T_critical, "Cold side temperature specification needs to be below critical temperature"

        h_condensate = PSI('H', 'Q', 0, 'T', T_hot + 273.15, working_fluid)
        p_condensate = PSI('P', 'Q', 0, 'T', T_hot + 273.15, working_fluid)
        h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        cd_va.set_attr(p=p_condensate / 1e5, h=h_condensate / 1e3)
        ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3)

    h_after_compression = h_evaporation + 50e3
    cp_cc.set_attr(h=h_after_compression / 1e3)

    nw.solve('design')

    return nw


def extract_plotting_results(nw):

    result_dict = {}
    result_dict.update(
        {cp.label: cp.get_plotting_data()[1] for cp in nw.comps['object']
         if cp.get_plotting_data() is not None})

    return result_dict


tespy_stable_nw = run_simple_heat_pump_model('CO2', "transcritical", 85, -20)
c1, c3, c4 = tespy_stable_nw.get_conn(['1', '3', '4'])
comp = tespy_stable_nw.get_comp('compressor')

# unset stable start value specifications and set target values instead
# for example compressor efficiency
c4.set_attr(h=None)
comp.set_attr(eta_s=0.9)

# or actual temperature levels instead of enthalpy
c1.set_attr(h=None, T=85)
# in transcritical case: set evaporation temperature and state of saturated
# gas after evaporation
c3.set_attr(h=None, p=None, x=1, T=-40)

tespy_stable_nw.solve('design')

tespy_stable_nw.print_results()

Just for curiosity reasons: Have you any specific plans with the software? How did you find it? :)

I'd be happy to hear from you in the future!

Have a nice day

matlaw · 2022-08-21T04:29:01Z

matlaw
Aug 21, 2022
Author

Thanks yet again, Francesco! Very, very helpful. I do have one more question. How would I change the code if I wanted the evaporation process to be isochoric instead of isobaric? I'm assuming I'd do something like calculate density after expansion, then use that as the input to the h_evaporation and p_evaporation calculations (instead of 'T' input). But I can't quite figure out how to make that change within the code. I'm not sure if there's something about the components as defined that makes this not possible. p_critical = PSI('p_critical', working_fluid) p_high = p_critical * 2 h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid) h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid) p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid) cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3) ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3) Matthew

…

On Tue, 16 Aug 2022 at 22:56, Francesco Witte ***@***.***> wrote: Sorry for the lengthy wait @matlaw <https://github.com/matlaw>, here you go: - For transcritical case you want to use the same specifications as supercritical for heat discharge and as subcritical for evaporation - When changing up parameters given the stable starting values, set fully saturated gas after evaporation and desired temperature level of evaporation instead of pressure and enthalpy from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Valve)from tespy.connections import Connectionfrom tespy.networks import Networkfrom CoolProp.CoolProp import PropsSI as PSI def run_simple_heat_pump_model(working_fluid, mode, T_hot, T_cold): nw = Network([working_fluid], T_unit='C', p_unit='bar', h_unit='kJ / kg') nw.set_attr(iterinfo=False) cp = Compressor('compressor') cc = CycleCloser('cycle_closer') cd = HeatExchangerSimple('condenser') va = Valve('expansion valve') ev = HeatExchangerSimple('evaporator') cc_cd = Connection(cc, 'out1', cd, 'in1', label='0') cd_va = Connection(cd, 'out1', va, 'in1', label='1') va_ev = Connection(va, 'out1', ev, 'in1', label='2') ev_cp = Connection(ev, 'out1', cp, 'in1', label='3') cp_cc = Connection(cp, 'out1', cc, 'in1', label='4') nw.add_conns(cc_cd, cd_va, va_ev, ev_cp, cp_cc) cd.set_attr(pr=1, Q=-1e6) ev.set_attr(pr=1) cc_cd.set_attr(fluid={working_fluid: 1}) if mode == "supercritical": # check if temperature values are actually supercritical T_critical = PSI('T_critical', working_fluid) - 273.15 assert T_hot > T_critical, "Hot side temperature specification needs to be above critical temperature" assert T_cold > T_critical, "Cold side temperature specification needs to be above critical temperature" # for a good guess of pressure levels p_critical = PSI('p_critical', working_fluid) p_high = p_critical * 2 # you do not want to be exactly at the critical point, here CoolProp ofte has difficulties solving # the equation of state p_low = p_critical * 1.1 # for a good guess of enthalpy values h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid) h_after_heat_input = PSI('H', 'P', p_low, 'T', T_cold + 273.15, working_fluid) cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3) ev_cp.set_attr(p=p_low / 1e5, h=h_after_heat_input / 1e3) elif mode == "transcritical": # check if temperature values are actually supercritical T_critical = PSI('T_critical', working_fluid) - 273.15 assert T_hot > T_critical, "Hot side temperature specification needs to be above critical temperature" assert T_cold < T_critical, "Cold side temperature specification needs to be below critical temperature" # for a good guess of pressure levels p_critical = PSI('p_critical', working_fluid) p_high = p_critical * 2 h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid) h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid) p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid) cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3) ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3) else: # check if temperature values are actually supercritical T_critical = PSI('T_critical', working_fluid) - 273.15 assert T_hot < T_critical, "Hot side temperature specification needs to be below critical temperature" assert T_cold < T_critical, "Cold side temperature specification needs to be below critical temperature" h_condensate = PSI('H', 'Q', 0, 'T', T_hot + 273.15, working_fluid) p_condensate = PSI('P', 'Q', 0, 'T', T_hot + 273.15, working_fluid) h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid) p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid) cd_va.set_attr(p=p_condensate / 1e5, h=h_condensate / 1e3) ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3) h_after_compression = h_evaporation + 50e3 cp_cc.set_attr(h=h_after_compression / 1e3) nw.solve('design') return nw def extract_plotting_results(nw): result_dict = {} result_dict.update( {cp.label: cp.get_plotting_data()[1] for cp in nw.comps['object'] if cp.get_plotting_data() is not None}) return result_dict tespy_stable_nw = run_simple_heat_pump_model('CO2', "transcritical", 85, -20)c1, c3, c4 = tespy_stable_nw.get_conn(['1', '3', '4'])comp = tespy_stable_nw.get_comp('compressor') # unset stable start value specifications and set target values instead# for example compressor efficiencyc4.set_attr(h=None)comp.set_attr(eta_s=0.9) # or actual temperature levels instead of enthalpyc1.set_attr(h=None, T=85)# in transcritical case: set evaporation temperature and state of saturated# gas after evaporationc3.set_attr(h=None, p=None, x=1, T=-40) tespy_stable_nw.solve('design') tespy_stable_nw.print_results() Just for curiosity reasons: Have you any specific plans with the software? How did you find it? :) I'd be happy to hear from you in the future! Have a nice day — Reply to this email directly, view it on GitHub <#352 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADGHYK5SXOK6QZV7DGQKIZTVZR5JZANCNFSM56PH4MEA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

2 replies

fwitte Aug 21, 2022
Maintainer

Dear Matthew,

But I can't quite figure out how to make that change within the code. I'm not sure if there's something about the components as defined that makes this not possible.

the heat exchanger implementation is directed to isobaric state of change (or with pressure losses). There would be three different ways, to make isochoric heat transfer work:

Using user defined equations API (a little bit tricky, but possible with current stable release):

from tespy.tools.fluid_properties import v_mix_ph
from tespy.tools.fluid_properties import dv_mix_dph
from tespy.tools.fluid_properties import dv_mix_pdh
from tespy.tools.helpers import UserDefinedEquation

def isochoric(self):
    return (
        v_mix_ph(self.conns[0].get_flow()) -
        v_mix_ph(self.conns[1].get_flow())
    )

def isochoric_deriv(self):
    # derivative to p for connection 0 in the list (c2 from list below)
    self.jacobian[self.conns[0]][1] = dv_mix_dph(self.conns[0].get_flow())
    # derivative to h for connection 0 in the list (c2 from list below)
    self.jacobian[self.conns[0]][2] = dv_mix_pdh(self.conns[0].get_flow())
    # derivative to p for connection 1 in the list (c3 from list below)
    self.jacobian[self.conns[1]][1] = -dv_mix_dph(self.conns[1].get_flow())
    # derivative to h for connection 1 in the list (c3 from list below)
    self.jacobian[self.conns[1]][2] = -dv_mix_pdh(self.conns[1].get_flow())
    return self.jacobian


ude = UserDefinedEquation(
    'isochoric evaporation', isochoric, isochoric_deriv, [c2, c3]
)
tespy_stable_nw.add_ude(ude)

eva = tespy_stable_nw.get_comp('evaporator')
# unset pressure ratio of evaporator
eva.set_attr(pr=None)

Using the Ref API (easiest, however not with current stable release since I only implemted it after your question ;)).

c2 = tespy_stable_nw.get_conn('2')
tespy_stable_nw.get_comp('evaporator').set_attr(pr=None)
c3.set_attr(v=Ref(c2, 1, 0))

Do you know how to install the developer version?

The components do have physical feasibility checks, since the "expected" behavior of the HeatExchangerSimple is isobaric heat transfer, a warning will be given after the simulation, since you have a slight increase of pressure. You can safely ignore the warning in this case, it might however be bothering to you. To disable it specifically for the evaporator, you can do the following:

# get rid of warnings
eva.get_attr("pr").max_val = 2
eva.get_attr("zeta").min_val = -1e12

Full script featuring the user defined equation API, since it is available on stable state of the software:

from tespy.components import (Compressor, CycleCloser, HeatExchangerSimple, Valve)
from tespy.connections import Connection, Ref
from tespy.networks import Network
from CoolProp.CoolProp import PropsSI as PSI


def run_simple_heat_pump_model(working_fluid, mode, T_hot, T_cold):
    nw = Network([working_fluid], T_unit='C', p_unit='bar', h_unit='kJ / kg')
    nw.set_attr(iterinfo=False)
    cp = Compressor('compressor')
    cc = CycleCloser('cycle_closer')
    cd = HeatExchangerSimple('condenser')
    va = Valve('expansion valve')
    ev = HeatExchangerSimple('evaporator')

    cc_cd = Connection(cc, 'out1', cd, 'in1', label='0')
    cd_va = Connection(cd, 'out1', va, 'in1', label='1')
    va_ev = Connection(va, 'out1', ev, 'in1', label='2')
    ev_cp = Connection(ev, 'out1', cp, 'in1', label='3')
    cp_cc = Connection(cp, 'out1', cc, 'in1', label='4')

    nw.add_conns(cc_cd, cd_va, va_ev, ev_cp, cp_cc)

    cd.set_attr(pr=1, Q=-1e6)
    ev.set_attr(pr=1)

    cc_cd.set_attr(fluid={working_fluid: 1})

    if mode == "supercritical":
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot > T_critical, "Hot side temperature specification needs to be above critical temperature"
        assert T_cold > T_critical, "Cold side temperature specification needs to be above critical temperature"

        # for a good guess of pressure levels
        p_critical = PSI('p_critical', working_fluid)
        p_high = p_critical * 2
        # you do not want to be exactly at the critical point, here CoolProp ofte has difficulties solving
        # the equation of state
        p_low = p_critical * 1.1

        # for a good guess of enthalpy values
        h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid)
        h_after_heat_input = PSI('H', 'P', p_low, 'T', T_cold + 273.15, working_fluid)

        cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3)
        ev_cp.set_attr(p=p_low / 1e5, h=h_after_heat_input / 1e3)
    elif mode == "transcritical":
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot > T_critical, "Hot side temperature specification needs to be above critical temperature"
        assert T_cold < T_critical, "Cold side temperature specification needs to be below critical temperature"

        # for a good guess of pressure levels
        p_critical = PSI('p_critical', working_fluid)
        p_high = p_critical * 2
        h_after_heat_discharge = PSI('H', 'P', p_high, 'T', T_hot + 273.15, working_fluid)
        h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        cd_va.set_attr(p=p_high / 1e5, h=h_after_heat_discharge / 1e3)
        ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3)

    else:
        # check if temperature values are actually supercritical
        T_critical = PSI('T_critical', working_fluid) - 273.15
        assert T_hot < T_critical, "Hot side temperature specification needs to be below critical temperature"
        assert T_cold < T_critical, "Cold side temperature specification needs to be below critical temperature"

        h_condensate = PSI('H', 'Q', 0, 'T', T_hot + 273.15, working_fluid)
        p_condensate = PSI('P', 'Q', 0, 'T', T_hot + 273.15, working_fluid)
        h_evaporation = PSI('H', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        p_evaporation = PSI('P', 'Q', 1, 'T', T_cold + 273.15, working_fluid)
        cd_va.set_attr(p=p_condensate / 1e5, h=h_condensate / 1e3)
        ev_cp.set_attr(p=p_evaporation / 1e5, h=h_evaporation / 1e3)

    h_after_compression = h_evaporation + 50e3
    cp_cc.set_attr(h=h_after_compression / 1e3)

    nw.solve('design')

    return nw


def extract_plotting_results(nw):

    result_dict = {}
    result_dict.update(
        {cp.label: cp.get_plotting_data()[1] for cp in nw.comps['object']
         if cp.get_plotting_data() is not None})

    return result_dict


tespy_stable_nw = run_simple_heat_pump_model('CO2', "transcritical", 85, -20)
c1, c2, c3, c4 = tespy_stable_nw.get_conn(['1', '2', '3', '4'])
comp = tespy_stable_nw.get_comp('compressor')

# unset stable start value specifications and set target values instead
# for example compressor efficiency
c4.set_attr(h=None)
comp.set_attr(eta_s=0.9)

# or actual temperature levels instead of enthalpy
c1.set_attr(h=None, T=85)
# in transcritical case: set evaporation temperature and state of saturated
# gas after evaporation
c3.set_attr(h=None, p=None, x=1, T=-40)

from tespy.tools.fluid_properties import v_mix_ph
from tespy.tools.fluid_properties import dv_mix_dph
from tespy.tools.fluid_properties import dv_mix_pdh
from tespy.tools.helpers import UserDefinedEquation

def isochoric(self):
    return (
        v_mix_ph(self.conns[0].get_flow()) -
        v_mix_ph(self.conns[1].get_flow())
    )

def isochoric_deriv(self):
    # derivative to p for connection 0 in the list (c2 from list below)
    self.jacobian[self.conns[0]][1] = dv_mix_dph(self.conns[0].get_flow())
    # derivative to h for connection 0 in the list (c2 from list below)
    self.jacobian[self.conns[0]][2] = dv_mix_pdh(self.conns[0].get_flow())
    # derivative to p for connection 1 in the list (c3 from list below)
    self.jacobian[self.conns[1]][1] = -dv_mix_dph(self.conns[1].get_flow())
    # derivative to h for connection 1 in the list (c3 from list below)
    self.jacobian[self.conns[1]][2] = -dv_mix_pdh(self.conns[1].get_flow())
    return self.jacobian


ude = UserDefinedEquation(
    'isochoric evaporation', isochoric, isochoric_deriv, [c2, c3]
)
tespy_stable_nw.add_ude(ude)

eva = tespy_stable_nw.get_comp('evaporator')
# unset pressure ratio of evaporator
eva.set_attr(pr=None)
# get rid of warnings
eva.get_attr("pr").max_val = 2
eva.get_attr("zeta").min_val = -1e12

tespy_stable_nw.solve('design')

tespy_stable_nw.print_results()

The user defined equation is actually a very potent functionality, you can learn more about it here: https://tespy.readthedocs.io/en/main/tespy_modules.html#user-defined-equations. It is possible to write down any kind of custom equation in your system as long as it is using the connection variables and you can correctly define the derivatives. You can also use the numeric_deriv method, if the derivatives cannot be written analytically or you just do not care and want quick solution ;).

fwitte Aug 21, 2022
Maintainer

To add: v_mix_ph calculates the specific volume of the mixture (or pure fluid). The connection property v, for example in c1.v.val, however referes to the volumetric flow, which would be the result of v_mix_ph(c.get_flow()) * c.m.val_SI. The actual specific volume in the results can be obtained from c1.v.val_SI / c1.m.val_SI or directly accessed with the vol keyword, so c1.vol.val.

While writing this down I find this very unintuitive, so I should change that sometime in the future.

matlaw · 2022-08-30T03:39:47Z

matlaw
Aug 30, 2022
Author

Apologies for not replying sooner - I'm very appreciative of your help. Unfortunately I haven't had a chance to look at it yet but I definitely will! This is sort of a side project at the moment, so it had to take the back seat for a little while. But I do really appreciate your help. Matthew

…

On Sun, 21 Aug 2022 at 03:50, Francesco Witte ***@***.***> wrote: To add: v_mix_ph calculates the specific volume of the mixture (or pure fluid). The connection property v, for example in c1.v.val, however referes to the volumetric flow, which would be the result of v_mix_ph(c.get_flow()) * c.m.val_SI. The actual specific volume in the results can be obtained from c1.v.val_SI / c1.m.val_SI or directly accessed with the vol keyword, so c1.vol.val. While writing this down I find this very unintuitive, so I should change that sometime in the future. — Reply to this email directly, view it on GitHub <#352 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADGHYK5ZITTEWL3IME3A2Q3V2ICYDANCNFSM56PH4MEA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

0 replies

matlaw · 2022-09-07T21:39:09Z

matlaw
Sep 7, 2022
Author

Hi Francesco, I implemented your code for the isochoric evaporation and also added in FluidDiagramProperty from fluprodia to map out the cycle paths on a T-s chart, and that all looks as I would expect. I've attached a copy of that code. However, I wrote some simple code to just check everything (see below) and I'm a bit confused. When I put in the temperature/pressures at the cycle start/endpoints given by the testpytranscritisochor results, I get the same enthalpy values as shown in the testpytranscritisochor results table: ##### RESULTS (Connection) ##### +----+-----------+-----------+-----------+------------+ | | m | p | h | T | |----+-----------+-----------+-----------+------------| | 0 | 5.281e+00 | 1.475e+02 | 5.556e+02 | 1.553e+02 | | 1 | 5.281e+00 | 1.475e+02 | 3.662e+02 | 6.500e+01 | | 2 | 5.281e+00 | 1.316e+01 | 3.662e+02 | -3.239e+01 | | 3 | 5.281e+00 | 1.683e+01 | 4.371e+02 | -2.500e+01 | | 4 | 5.281e+00 | 1.475e+02 | 5.556e+02 | 1.553e+02 | +----+-----------+-----------+-----------+------------+ But it looks to me like testpytranscritisochor.py is using the change in enthalpy over the isochoric process when calculating evaporator Qin, but shouldn't it be using change in internal energy since the pressure remains fixed? Here is the heat exchange part of the testpytrancritisochor output: ##### RESULTS (HeatExchangerSimple) ##### +------------+-----------+----------+-----------+-----+-----+------+------+--------+ | | Q | pr | zeta | D | L | ks | kA | Tamb | |------------+-----------+----------+-----------+-----+-----+------+------+--------| | condenser | -1.00e+06 | 1.00e+00 | 0.00e+00 | nan | nan | nan | nan | nan | | evaporator | 3.74e+05 | 1.28e+00 | -7.11e+05 | nan | nan | nan | nan | nan | +------------+-----------+----------+-----------+-----+-----+------+------+--------+ When I calculate from evaporator Q from deltaU I get Qin = 329.6kJ (vs 373.7kJ if use deltaH). I'm just wondering if something in the simple heat exchanger model needs to be adjusted so that delta internal energy is used when calculating evaporator Q during isochoric? Or am I mistaken about how to calculate Qin during isochoric evaporation? When I change both models to isobaric evaporation, I get matching results. Here is my simple code that I use to check results: from CoolProp.CoolProp import PropsSI import numpy as np from prettytable import PrettyTable import matplotlib.pyplot as plt #Assume temperature and quality prior to compression (right at saturated vapor line?) T1 = -25 + 273.15 #temperature in K Q1 = 1 #quality #calculate pressure, enthalpy from fluid properties P1 = PropsSI('P','T',T1,'Q',Q1,'CO2') #Pressure in Pa H1 = PropsSI('H','T',T1,'Q',Q1,'CO2') #Enthalpy in J/kg rho1 = PropsSI('D','T',T1,'Q',Q1,'CO2') #density in kg/m3 U1 = PropsSI('U','T',T1,'Q',Q1,'CO2') #internal energy in J/kg #Assume compressor work and CO2 mass flow are known Wcomp = 626.058E3 #work in watts mdot = 5.281 #C02 mass flow in kg/s #calculate enthalpy after compression H2 = Wcomp / mdot + H1 #assume fluid pressure after compression is known P2 = 147.5 * 100 * 1000 #pressure in Pa #calculate T2 from fluid propertioes T2 = PropsSI('T','P',P2,'H',H2,'CO2') U2 = PropsSI('U','P',P2,'H',H2,'CO2') Q2 = PropsSI('Q','P',P2,'H',H2,'CO2') rho2 = PropsSI('D', 'P', P2, 'H', H2, 'CO2') #assume fluid cooling is isobaric P3 = P2 #start with a guess for T3 T3 = 65 + 273.15 #calculate H3 from fluid properties H3=PropsSI('H','P',P3,'T',T3,'CO2') #calculate Qout during fluid cooling Qout=mdot * (H3 - H2) #calculate enthalpy after expansion H4=H3 #evaporation is isochoric process - can add a delta term if required rho4=rho1 #calculate internal energy from fluid properties U4=PropsSI('U','H',H4,'D',rho4,'CO2') #internal energy in J/kg #calculate Qin during evaporation Qin=mdot * (U1 - U4) Thanks again for your help! Matthew

…

On Mon, 29 Aug 2022 at 20:39, Matthew Lawrence ***@***.***> wrote: Apologies for not replying sooner - I'm very appreciative of your help. Unfortunately I haven't had a chance to look at it yet but I definitely will! This is sort of a side project at the moment, so it had to take the back seat for a little while. But I do really appreciate your help. Matthew On Sun, 21 Aug 2022 at 03:50, Francesco Witte ***@***.***> wrote: > To add: v_mix_ph calculates the specific volume of the mixture (or pure > fluid). The connection property v, for example in c1.v.val, however > referes to the volumetric flow, which would be the result of > v_mix_ph(c.get_flow()) * c.m.val_SI. The actual specific volume in the > results can be obtained from c1.v.val_SI / c1.m.val_SI or directly > accessed with the vol keyword, so c1.vol.val. > > While writing this down I find this very unintuitive, so I should change > that sometime in the future. > > — > Reply to this email directly, view it on GitHub > <#352 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/ADGHYK5ZITTEWL3IME3A2Q3V2ICYDANCNFSM56PH4MEA> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >

0 replies

matlaw · 2022-09-08T22:34:39Z

matlaw
Sep 8, 2022
Author

Sorry - I just realized a couple of things: 1) I sent you the wrong python code. I've attached the right one 2) I think I figured out my confusion. Since the evaporation process is isochoric, part of the energy will go to maintaining constant volume. So the Q calculated by your code would be correct (and gives energy balance for the cycle), but the cooling applied to the environment outside the heat exchanger would be less (and given by mdot deltaU instead of mdot deltaH). I think. Matthew

…

On Wed, 7 Sept 2022 at 14:38, Matthew Lawrence ***@***.***> wrote: Hi Francesco, I implemented your code for the isochoric evaporation and also added in FluidDiagramProperty from fluprodia to map out the cycle paths on a T-s chart, and that all looks as I would expect. I've attached a copy of that code. However, I wrote some simple code to just check everything (see below) and I'm a bit confused. When I put in the temperature/pressures at the cycle start/endpoints given by the testpytranscritisochor results, I get the same enthalpy values as shown in the testpytranscritisochor results table: ##### RESULTS (Connection) ##### +----+-----------+-----------+-----------+------------+ | | m | p | h | T | |----+-----------+-----------+-----------+------------| | 0 | 5.281e+00 | 1.475e+02 | 5.556e+02 | 1.553e+02 | | 1 | 5.281e+00 | 1.475e+02 | 3.662e+02 | 6.500e+01 | | 2 | 5.281e+00 | 1.316e+01 | 3.662e+02 | -3.239e+01 | | 3 | 5.281e+00 | 1.683e+01 | 4.371e+02 | -2.500e+01 | | 4 | 5.281e+00 | 1.475e+02 | 5.556e+02 | 1.553e+02 | +----+-----------+-----------+-----------+------------+ But it looks to me like testpytranscritisochor.py is using the change in enthalpy over the isochoric process when calculating evaporator Qin, but shouldn't it be using change in internal energy since the pressure remains fixed? Here is the heat exchange part of the testpytrancritisochor output: ##### RESULTS (HeatExchangerSimple) ##### +------------+-----------+----------+-----------+-----+-----+------+------+--------+ | | Q | pr | zeta | D | L | ks | kA | Tamb | |------------+-----------+----------+-----------+-----+-----+------+------+--------| | condenser | -1.00e+06 | 1.00e+00 | 0.00e+00 | nan | nan | nan | nan | nan | | evaporator | 3.74e+05 | 1.28e+00 | -7.11e+05 | nan | nan | nan | nan | nan | +------------+-----------+----------+-----------+-----+-----+------+------+--------+ When I calculate from evaporator Q from deltaU I get Qin = 329.6kJ (vs 373.7kJ if use deltaH). I'm just wondering if something in the simple heat exchanger model needs to be adjusted so that delta internal energy is used when calculating evaporator Q during isochoric? Or am I mistaken about how to calculate Qin during isochoric evaporation? When I change both models to isobaric evaporation, I get matching results. Here is my simple code that I use to check results: from CoolProp.CoolProp import PropsSI import numpy as np from prettytable import PrettyTable import matplotlib.pyplot as plt #Assume temperature and quality prior to compression (right at saturated vapor line?) T1 = -25 + 273.15 #temperature in K Q1 = 1 #quality #calculate pressure, enthalpy from fluid properties P1 = PropsSI('P','T',T1,'Q',Q1,'CO2') #Pressure in Pa H1 = PropsSI('H','T',T1,'Q',Q1,'CO2') #Enthalpy in J/kg rho1 = PropsSI('D','T',T1,'Q',Q1,'CO2') #density in kg/m3 U1 = PropsSI('U','T',T1,'Q',Q1,'CO2') #internal energy in J/kg #Assume compressor work and CO2 mass flow are known Wcomp = 626.058E3 #work in watts mdot = 5.281 #C02 mass flow in kg/s #calculate enthalpy after compression H2 = Wcomp / mdot + H1 #assume fluid pressure after compression is known P2 = 147.5 * 100 * 1000 #pressure in Pa #calculate T2 from fluid propertioes T2 = PropsSI('T','P',P2,'H',H2,'CO2') U2 = PropsSI('U','P',P2,'H',H2,'CO2') Q2 = PropsSI('Q','P',P2,'H',H2,'CO2') rho2 = PropsSI('D', 'P', P2, 'H', H2, 'CO2') #assume fluid cooling is isobaric P3 = P2 #start with a guess for T3 T3 = 65 + 273.15 #calculate H3 from fluid properties H3=PropsSI('H','P',P3,'T',T3,'CO2') #calculate Qout during fluid cooling Qout=mdot * (H3 - H2) #calculate enthalpy after expansion H4=H3 #evaporation is isochoric process - can add a delta term if required rho4=rho1 #calculate internal energy from fluid properties U4=PropsSI('U','H',H4,'D',rho4,'CO2') #internal energy in J/kg #calculate Qin during evaporation Qin=mdot * (U1 - U4) Thanks again for your help! Matthew On Mon, 29 Aug 2022 at 20:39, Matthew Lawrence ***@***.***> wrote: > Apologies for not replying sooner - I'm very appreciative of your help. > Unfortunately I haven't had a chance to look at it yet but I > definitely will! > > This is sort of a side project at the moment, so it had to take the back > seat for a little while. But I do really appreciate your help. > > Matthew > > On Sun, 21 Aug 2022 at 03:50, Francesco Witte ***@***.***> > wrote: > >> To add: v_mix_ph calculates the specific volume of the mixture (or pure >> fluid). The connection property v, for example in c1.v.val, however >> referes to the volumetric flow, which would be the result of >> v_mix_ph(c.get_flow()) * c.m.val_SI. The actual specific volume in the >> results can be obtained from c1.v.val_SI / c1.m.val_SI or directly >> accessed with the vol keyword, so c1.vol.val. >> >> While writing this down I find this very unintuitive, so I should change >> that sometime in the future. >> >> — >> Reply to this email directly, view it on GitHub >> <#352 (reply in thread)>, >> or unsubscribe >> <https://github.com/notifications/unsubscribe-auth/ADGHYK5ZITTEWL3IME3A2Q3V2ICYDANCNFSM56PH4MEA> >> . >> You are receiving this because you were mentioned.Message ID: >> ***@***.***> >> >

1 reply

fwitte Sep 19, 2022
Maintainer

@matlaw, I am sorry for the late reply. I am on holidays currently, will come back to this soon. Best

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

no simple heat pump model solution when changing fluid from NH3 to CO2 #352

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 5 comments 6 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

no simple heat pump model solution when changing fluid from NH3 to CO2 #352

matlaw Aug 14, 2022

Replies: 5 comments · 6 replies

fwitte Aug 14, 2022 Maintainer

matlaw Aug 15, 2022 Author

fwitte Aug 15, 2022 Maintainer

fwitte Aug 17, 2022 Maintainer

matlaw Aug 21, 2022 Author

fwitte Aug 21, 2022 Maintainer

fwitte Aug 21, 2022 Maintainer

matlaw Aug 30, 2022 Author

matlaw Sep 7, 2022 Author

matlaw Sep 8, 2022 Author

fwitte Sep 19, 2022 Maintainer

matlaw
Aug 14, 2022

Replies: 5 comments 6 replies

fwitte
Aug 14, 2022
Maintainer

matlaw Aug 15, 2022
Author

fwitte Aug 15, 2022
Maintainer

fwitte Aug 17, 2022
Maintainer

matlaw
Aug 21, 2022
Author

fwitte Aug 21, 2022
Maintainer

fwitte Aug 21, 2022
Maintainer

matlaw
Aug 30, 2022
Author

matlaw
Sep 7, 2022
Author

matlaw
Sep 8, 2022
Author

fwitte Sep 19, 2022
Maintainer