Heat pump issue (Cooling mode) #528

rafaelmarqferreira · 2024-07-19T08:46:05Z

rafaelmarqferreira
Jul 19, 2024

Hey there. I´ve been working on a heat pump model to test the feasibility of a project which would lead to building a test bench later on. The test bench would have a 4 ways reversing valve, which would enable the possibility to work in both heating and cooling modes. While modelling the heating mode, everything went smoothly, but when it came to cooling mode, linear dependencies errors seem to show up, even though I don´t believe to be over defining the model. I´ve left it down below if need.

Is it me overseeing something here?

Thanks in advance!

from tespy.networks import Network
from CoolProp.CoolProp import PropsSI as PSI
import matplotlib.pyplot as plt
from fluprodia import FluidPropertyDiagram
import math
from scipy import constants
import sys
from scipy.special import iv, kv
from tespy.connections import Connection
from tespy.components import (CycleCloser, Compressor, Valve, Pump, HeatExchanger, Source, Sink)

heatPump_heating = Network(T_unit='C', p_unit='Pa', h_unit='kJ / kg', iterinfo=False)

closer = CycleCloser('cycle closer') #Chamar componente ciclo de fecho
Cooling1 = HeatExchanger('Cooling 1')
Cooling2 = HeatExchanger('Cooling 2') #Chamar componente permutador de calor
valve = Valve('expansion valve') #Chamar componente válvula
compressor = Compressor('compressor') #Chamar componente compressor
pump = Pump('pump') #Chamar componente bomba

source1 = Source('ambIn1')
source2 = Source('ambIn2')
sink1 = Sink('ambOut1')
sink2 = Sink('ambOut2')

c0 = Connection(Cooling2, 'out2', compressor, 'in1', label= 'Estado inicial')
c1 = Connection(compressor, 'out1', Cooling1, 'in1', label= 'Pós compressão')
c2 = Connection(Cooling1, "out1", closer, "in1", label= "Pós arrefecimento")
c3 = Connection(closer, "out1", valve, "in1", label= "Cycle Closer")
c4 = Connection(valve, "out1", Cooling2, "in2", label= "Pós expansão")

c10 = Connection(source1, "out1", Cooling1, "in2", label= "Ambiente")
c11 = Connection(Cooling1, "out2", sink1, "in1", label= "Ambiente aquecido")

c20 = Connection(source2, "out1", Cooling2, "in1", label= "Água quente")
c21 = Connection(Cooling2, "out1", sink2, "in1", label= "Água fria")


heatPump_heating.add_conns(c0, c1, c2, c3, c4, c10, c11, c20, c21)

wf = "R32"
massRateWF = 0.037 #kg/s
T_wf_in = -33
pressureIn = 1.8e5
pressureOut = 16e5
h0 = PSI("H", "P", pressureIn, "T", 273.15 + T_wf_in, wf) / 1000
fluidCheck = PSI("T", "P", pressureIn, "Q", 1, wf) - 273.15
T_satWF = PSI("T", "P", pressureOut,"Q", 1 , wf) - 273.15
T_amb = -28.3

T_cooling_in2 = 8 #Água
coolingHeated1 = 25 #Ar
coolingHeated2 = 5 #Água
coolingVolumeRate2 = 2 / 3600 #m3/s
coolingFluid2 = "Water"
coolingFluid1 = "Air"
coolingDensity2 = PSI("D", "T", 273.15 + T_cooling_in2, "P", 3e5, coolingFluid2) #kg/m3 Água
massRateCooling2 = coolingVolumeRate2 * coolingDensity2 #kg/s
coolingVolume1 = 800 / 3600 #m3/s
massRateCooling1 = PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * coolingVolume1

t_pack = T_amb
t_pack_final = 20
m_pack = 335
cp_pack = 500
a_pack = 1.073 #m
b_pack = 0.578 #m
c_pack = 0.586 #m
volume_pack = a_pack * b_pack * c_pack #m3
transfer_area_pack = 4 * a_pack * c_pack + a_pack * b_pack #m2

if T_wf_in <= fluidCheck:
    print("Existência de líquido / mistura bifásica no compressor. Aumente a temperatura de entrada.")
else:
    print("Modo de aquecimento")
    condenserPressureLossRatioHot = 1
    evaporatorPressureLossRatioHot = 1

    efficiencyCompressor = 0.85

    efficiencyPump = 0.65

    Cooling1.set_attr(pr1 = condenserPressureLossRatioHot, pr2 = 1)
    compressor.set_attr(eta_s = efficiencyCompressor, pr = pressureOut / pressureIn)
    Cooling2.set_attr(pr1 = evaporatorPressureLossRatioHot)

    valve.set_attr(pr = pressureIn / pressureOut)

    c0.set_attr(fluid={wf: 1}, m = massRateWF, T = T_wf_in, p = pressureIn)

    c10.set_attr(fluid={coolingFluid1: 1}, p = 1.01325e5, T = T_amb, m = massRateCooling1)
    c11.set_attr(T = coolingHeated1)

    c20.set_attr(fluid={coolingFluid2: 1}, p = 3e5, T = T_cooling_in2)
    c21.set_attr(T = coolingHeated2)
    
    heatPump_heating.solve(mode='design')

    heatPump_heating.save('hp_design_heating')
    heatPump_heating.print_results()

    cop = abs(Cooling1.Q.val) / compressor.P.val

    result_dict = {}
    result_dict.update({Cooling2.label: Cooling2.get_plotting_data()[2]})
    result_dict.update({compressor.label: compressor.get_plotting_data()[1]})
    result_dict.update({Cooling1.label: Cooling1.get_plotting_data()[1]})
    result_dict.update({valve.label: valve.get_plotting_data()[1]})

    diagram = FluidPropertyDiagram(wf)
    diagram.set_unit_system(T='°C', p='Pa', h='kJ/kg')

    for key, data in result_dict.items():
        result_dict[key]['datapoints'] = diagram.calc_individual_isoline(**data)

    diagram.calc_isolines()

    fig, ax = plt.subplots(1, figsize=(16, 10))
    diagram.draw_isolines(fig, ax, 'logph', x_min= c3.h.val - 350, x_max= c1.h.val + 100, y_min= 1e4  , y_max= c1.p.val * 10)

    for key in result_dict.keys():
        datapoints = result_dict[key]['datapoints']
        ax.plot(datapoints['h'], datapoints['p'], color='#ff0000')
        ax.scatter(datapoints['h'][0], datapoints['p'][0], color='#ff0000')

    print('O COP da bomba de calor com ' + str(wf)+ ' é ' + str(cop))
    fig.savefig(str(wf) + '_logph_heating.svg')

    #Dimensionamento do permutador de placas paralelas

    port_D = 0.016 #m
    plate_thickness = 0.0005 #m / d_p
    plate_conductivity = 14.9 # W/m.K
    chevron_angle = 60 # º
    l_v = 0.250 #m
    l_p = l_v - port_D #m
    w_p = 0.111 #m
    l_h = 0.050 #m
    d_g = 0.00236 # Espaço médio entre placas / channel gap / b #m
    # corrugation_pitch = 0.0075 #
    plates = 50
    channels_refrigerant = (plates - 2) / 2
    channels_cooling = plates / 2
    area_catalogo = 0.0026 * (plates - 2) #m2
    compressed_pitch = d_g + plate_thickness #m
    amplitude = d_g / 2 #m
    wave_length = 0.008 #m experimentação
    omega = (math.pi * d_g) / wave_length #m
    enlargement_factor = (1 / 6) * (1 + math.sqrt(1 + (omega ** 2)) + 4 * math.sqrt(1 + ((omega ** 2) / 2)))
    d_h = 2 * d_g / enlargement_factor #m
    d_e = 2 * d_g #m
    corrugation_aspect_ratio = 2 * d_g / wave_length #m
    area_cross = d_g * w_p #m2
    

    if c4.x.val == -1:
        q1 = massRateWF * (PSI("H", "Q", 0, "P", pressureIn, wf) - c4.h.val * 1000) #W
        T_cooling_mid1 = T_cooling_in2 + q1 / (massRateCooling2 * PSI("C", "T", T_cooling_in2 + 273.15, "P", c20.p.val, coolingFluid2))
        c_h_1 = massRateWF * PSI("C", "T", 273.15 + ((PSI("T", "Q", 0, "P", pressureIn, wf) - 273.15) + c4.T.val) / 2, "P", pressureIn, wf) #W
        c_c_1 = massRateCooling2 * PSI("C", "T", 273.15 + (T_cooling_mid1 + c20.T.val) / 2, "P", c20.p.val, coolingFluid2) #W
        c_min_1 = min(c_h_1, c_c_1)
        c_max_1 = max(c_h_1, c_c_1)
        Q_1_max = c_min_1 * (c20.T.val - c4.T.val)
        efficiency_1 = q1 / Q_1_max
        Cr_1 = c_min_1 / c_max_1
        ntu_1 = (1 / (Cr_1 - 1)) * math.log((efficiency_1 - 1) / (efficiency_1 * Cr_1 - 1))

        q2 = massRateWF * (PSI("H", "Q", 1, "P", pressureIn, wf) - PSI("H", "Q", 0, "P", pressureIn, wf))
        T_cooling_mid2 = T_cooling_mid1 + q2 / (massRateCooling2 * PSI("C", "T", T_cooling_mid1 + 273.15, "P", c20.p.val, coolingFluid2))
        c_c_2 = massRateCooling2 * PSI("C", "T", 273.15 + (T_cooling_mid1 + T_cooling_mid2) / 2, "P", c20.p.val, coolingFluid2) #W
        Q_2_max = c_c_2 * (T_cooling_mid1 - (PSI("T", "Q", 0, "P", pressureIn, wf) - 273.15))
        efficiency_2 = q2 / Q_2_max
        ntu_2 = - math.log(1 - efficiency_2)
    else:
        q1 = 0
        ntu_1 = 0
        c_min_1 = 0
        q2 = massRateWF * (PSI("H", "Q", 1, "P", pressureIn, wf) - PSI("H", "Q", 0, "P", pressureIn, wf)) * (1 - c4.x.val)
        T_cooling_mid2 = T_cooling_in2
        T_cooling_mid1 = T_cooling_in2
        c_c_2 = massRateCooling2 * PSI("C", "T", 273.15 + T_cooling_mid2, "P", c20.p.val, coolingFluid2) #W
        Q_2_max = c_c_2 * (T_cooling_mid2 - (PSI("T", "Q", 0, "P", pressureIn, wf) - 273.15))
        efficiency_2 = q2 / Q_2_max
        ntu_2 = - math.log(1 - efficiency_2)

    q3 = abs(Cooling2.Q.val) - q2 - q1
    c_h_3 = massRateWF * PSI("C", "T", 273.15 + ((PSI("T", "Q", 1, "P", pressureIn, wf) - 273.15) + c0.T.val) / 2, "P", pressureIn, wf) #W
    c_c_3 = massRateCooling2 * PSI("C", "T", 273.15 + (T_cooling_mid2 + c21.T.val) / 2, "P", c20.p.val, coolingFluid2) #W
    c_min_3 = min(c_h_3, c_c_3)
    c_max_3 = max(c_h_3, c_c_3)
    Cr_3 = c_min_3 / c_max_3
    Q_3_max = c_min_3 * (T_cooling_mid2 - (PSI("T", "Q", 1, "P", pressureIn, wf) - 273.15))
    efficiency_3 = q3 / Q_3_max
    ntu_3 = (1 / (Cr_3 - 1)) * math.log((efficiency_3 - 1) / (efficiency_3 * Cr_3 - 1))

    G_port_h = 4 * c1.m.val / (math.pi * port_D)
    G_port_c = 4 * c20.m.val / (math.pi * port_D)

    A_flow_h = w_p * d_g * channels_refrigerant
    A_flow_c = w_p * d_g * channels_cooling

    G_channel_h = massRateWF / A_flow_h
    G_channel_c = massRateCooling2 / A_flow_c
    Re_channel_h1 = G_channel_h * d_h / PSI("V", "T", 273.15 + (c4.T.val + PSI("T", "P", pressureIn, "Q", 0, wf)) / 2, "P", pressureIn, wf)
    Re_channel_h3 = G_channel_h * d_h / PSI("V", "T", 273.15 + (PSI("T", "P", pressureIn, "Q", 0, wf) + T_wf_in) / 2, "P", pressureIn, wf)
    Re_channel_c1 = G_channel_c * d_h / PSI("V", "T", 273.15 + (c20.T.val + T_cooling_mid1) / 2, "P", c20.p.val, coolingFluid2)
    Re_channel_c2 = G_channel_c * d_h / PSI("V", "T", 273.15 + (T_cooling_mid1 + T_cooling_mid2) / 2, "P", c20.p.val, coolingFluid2)
    Re_channel_c3 = G_channel_c * d_h / PSI("V", "T", 273.15 + (T_cooling_mid2 + c21.T.val) / 2, "P", c20.p.val, coolingFluid2)

    R_f_h = 0.0001
    R_f_c = 0.0001

    if c4.x.val == -1:

        Pr_c1 = PSI("V", "T", 273.15 + (c20.T.val + c21.T.val) / 2, "P", c20.p.val, coolingFluid2) * PSI("C", "T", 273.15 + (c20.T.val + c21.T.val) / 2, "P", c20.p.val, coolingFluid2) / PSI("L", "T", 273.15 + (c20.T.val + c21.T.val) / 2, "P", c20.p.val, coolingFluid2)
        Pr_h1 = PSI("V", "T", 273.15 + (c4.T.val + PSI("T", "P", pressureIn, "Q", 0, wf)) / 2, "P", pressureIn, wf) * PSI("C", "T", 273.15 + (c4.T.val + PSI("T", "P", pressureIn, "Q", 0, wf)) / 2, "P", pressureIn, wf) / PSI("L", "T", 273.15 + (c4.T.val + PSI("T", "P", pressureIn, "Q", 0, wf)) / 2, "P", pressureIn, wf)
        h_c1 = 0.724 * (PSI("L", "T", 273.15 + (c20.T.val + c21.T.val) / 2, "P", c20.p.val, coolingFluid2) / d_h) * ((chevron_angle / 30) ** (0.646)) * (Re_channel_c1 ** 0.583) * (Pr_c1 ** (1/3))
        h_1 = 0.724 * (PSI("L", "T", 273.15 + (c4.T.val + PSI("T", "P", pressureIn, "Q", 0, wf)) / 2, "P", pressureIn, wf) / d_h) * ((chevron_angle / 30) ** (0.646)) * (Re_channel_h1 ** 0.583) * (Pr_h1 ** (1/3))
        u_1 = 1 / ((1 / h_c1) + (1 / h_1))
        a_1 = ntu_1 * c_min_1 / u_1
    else:
        h_1 = 0
        ntu_1 = 0
        a_1 = 0
    
    if c4.x.val == -1:
        average_quality_in_out = 0.5
    else:
        average_quality_in_out = (1 + c4.x.val) / 2

    We_m = (G_channel_h ** 2) * d_h / ((average_quality_in_out * PSI("D", "Q", 1, "P", pressureIn, wf) + (1 - average_quality_in_out) * PSI("D", "Q", 0, "P", pressureIn, wf)) * PSI("I", "Q", 0, "P", pressureIn, wf))
    Re_lo = G_channel_h * d_h / PSI("V", "Q", 0, "P", pressureIn, wf)
    Re_vo = G_channel_h * d_h / PSI("V", "Q", 1, "P", pressureIn, wf)
    Pho_ratio = PSI("D", "Q", 0, "P", pressureIn, wf) / PSI("D", "Q", 1, "P", pressureIn, wf)
    chevron_ratio = chevron_angle / 30
    Re_l = G_channel_h * (1 - average_quality_in_out) * d_h / PSI("V", "Q", 0, "P", pressureIn, wf)
    Re_v = G_channel_h * average_quality_in_out * d_h / PSI("V", "Q", 0, "P", pressureIn, wf)
    Bo = massRateWF * (PSI("H", "P", pressureIn, "Q", 1, wf) -  c4.h.val * 1000) / (area_catalogo * G_channel_h * (PSI("H", "P", pressureIn, "Q", 1, wf) - PSI("H", "P", pressureIn, "Q", 0, wf)))
    Bd = (PSI("D", "Q", 0, "P", pressureIn, wf) - PSI("D", "Q", 1, "P", pressureIn, wf)) * constants.g * (d_h ** 2) / PSI("I", "Q", 0, "P", pressureIn, wf)

    d_0 = 0.0146 * 35 * ((2 * PSI("I", "Q", 0, "P", pressureIn, wf) / (constants.g * (PSI("D", "Q", 0, "P", pressureIn, wf) - PSI("D", "Q", 1, "P", pressureIn, wf)))) ** 0.5)
    thermal_diffusivity = (PSI("L", "P", pressureIn, "T", 273.15 + 0.0000001 + c0.T.val, wf)) / (PSI("D", "Q", 0, "P", pressureIn, wf) * PSI("C", "Q", 0, "P", pressureIn, wf))
    Pr_l = PSI("D", "Q", 0, "P", pressureIn, wf) * PSI("V", "Q", 0, "P", pressureIn, wf) / PSI("L", "Q", 0, "P", pressureIn, wf)
    Nu_tp = 0.00187 * ((massRateWF * (PSI("H", "P", pressureIn, "Q", 1, wf) -  c4.h.val * 1000) * d_0 / (PSI("L", "P", pressureIn, "T", 273.15 + 0.0000001 + c0.T.val, wf))) ** 0.56) * (((PSI("H", "P", pressureIn, "Q", 1, wf) - PSI("H", "P", pressureIn, "Q", 0, wf)) * d_0 / (thermal_diffusivity ** 2)) ** 0.31) * (Pr_l ** 0.33)
    h_tp = Nu_tp * (PSI("L", "T", 273.15 + 0.001 + c0.T.val, "P", pressureIn, wf)) / l_p

    Pr_c2 = PSI("V", "T", 273.15 + (T_cooling_mid1 + T_cooling_mid2) / 2, "P", c20.p.val, coolingFluid2) * PSI("C", "T", 273.15 + (T_cooling_mid1 + T_cooling_mid2) / 2, "P", c20.p.val, coolingFluid2) / PSI("L", "T", 273.15 + (T_cooling_mid1 + T_cooling_mid2) / 2, "P", c20.p.val, coolingFluid2)
    h_c2 = 0.724 * (PSI("L", "T", 273.15 + (T_cooling_mid1 + T_cooling_mid2) / 2, "P", c20.p.val, coolingFluid2) / d_h) * ((chevron_angle / 30) ** (0.646)) * (Re_channel_c2 ** 0.583) * (Pr_c2 ** (1/3))
    u_2 = 1 / ((1 / h_c2) + (1 / h_tp))
    a_2 = ntu_2 * c_c_2 / u_2

    Pr_h3 = PSI("V", "T", 273.15 + (PSI("T", "P", pressureIn, "Q", 0, wf) + T_wf_in) / 2, "P", pressureIn, wf) * PSI("C", "T", 273.15 + (PSI("T", "P", pressureIn, "Q", 0, wf) + T_wf_in) / 2, "P", pressureIn, wf) / PSI("L", "T", 273.15 + (PSI("T", "P", pressureIn, "Q", 0, wf) + T_wf_in) / 2, "P", pressureIn, wf)
    h_3 = 0.724 * (PSI("L", "T", 273.15 + (PSI("T", "P", pressureIn, "Q", 0, wf) + T_wf_in) / 2, "P", pressureIn, wf) / d_h) * ((chevron_angle / 30) ** (0.646)) * (Re_channel_h3 ** 0.583) * (Pr_h3 ** (1/3))
    Pr_c3 = PSI("V", "T", 273.15 + (T_cooling_mid2 + c21.T.val) / 2, "P", c20.p.val, coolingFluid2) * PSI("C", "T", 273.15 + (T_cooling_mid2 + c21.T.val) / 2, "P", c20.p.val, coolingFluid2) / PSI("L", "T", 273.15 + (T_cooling_mid2 + c21.T.val) / 2, "P", c20.p.val, coolingFluid2)
    h_c3 = 0.724 * (PSI("L", "T", 273.15 + (T_cooling_mid2 + c21.T.val) / 2, "P", c20.p.val, coolingFluid2) / d_h) * ((chevron_angle / 30) ** (0.646)) * (Re_channel_c3 ** 0.583) * (Pr_c3 ** (1/3))
    u_3 = 1 / ((1 / h_c3) + (1 / h_3))
    a_3 = ntu_3 * c_min_3 / u_3

    area_needed = a_1 + a_2 + a_3
    
    F_r_f = 0.183 * (chevron_ratio ** 2) - 0.275 * chevron_ratio + 1.1
    pho_average = 1 / ((average_quality_in_out / PSI("D", "Q", 1, "P", pressureIn, wf)) + (1 - average_quality_in_out) / PSI("D", "Q", 0, "P", pressureIn, wf))
    v_tp = pho_average * (average_quality_in_out * PSI("V", "Q", 1, "P", pressureIn, wf) / PSI("D", "Q", 1, "P", pressureIn, wf) + (1 - average_quality_in_out) * PSI("V", "Q", 0, "P", pressureIn, wf) / PSI("D", "Q", 0, "P", pressureIn, wf))
    Re_tp = G_channel_h * d_h / v_tp

    f_c1 = 0.8 * (Re_channel_c1 ** (-0.25))
    f_c2 = 0.8 * (Re_channel_c2 ** (-0.25))
    f_c3 = 0.8 * (Re_channel_c3 ** (-0.25))
    f_h1 = 0.8 * (Re_channel_h1 ** (-0.25))
    f_tp = 38100 * F_r_f / ((Re_tp ** 0.9) * (Pho_ratio ** 0.16))
    f_h3 = 0.8 * (Re_channel_h3 ** (-0.25))

    delta_p_f_c1 = (4 * f_c1 * l_p * (G_channel_c ** 2)) / (2 * PSI("D", "P", c20.p.val, "T", 273.15 + T_cooling_mid1, coolingFluid2) * d_h)
    delta_p_f_c2 = (4 * f_c2 * l_p * (G_channel_c ** 2)) / (2 * PSI("D", "P", c20.p.val, "T", 273.15 + T_cooling_mid2, coolingFluid2) * d_h)
    delta_p_f_c3 = (4 * f_c3 * l_p * (G_channel_c ** 2)) / (2 * PSI("D", "P", c20.p.val, "T", 273.15 + (c21.T.val + T_cooling_mid2) / 2, coolingFluid2) * d_h)

    if c4.x.val == -1:
        delta_p_f_h1 = (4 * f_h1 * l_p * (G_channel_h ** 2)) / (2 * PSI("D", "P", pressureIn, "T", 273.15 + (c4.T.val + PSI("T", "P", pressureIn, "Q", 0, wf)) / 2, wf) * d_h)
    else:
        delta_p_f_h1 = 0
    
    delta_p_f_h2 = (4 * f_tp * l_p * (G_channel_h ** 2)) / ((PSI("D", "P", pressureIn, "Q", 1, wf) + PSI("D", "P", pressureIn, "Q", 0, wf)) * d_h)
    delta_p_f_h3 = (4 * f_h3 * l_p * (G_channel_h ** 2)) / (2 * PSI("D", "P", pressureIn, "T", 273.15 + (PSI("T", "P", pressureIn, "Q", 0, wf) + T_wf_in) / 2, wf) * d_h)

    delta_p_f_h = delta_p_f_h1 + delta_p_f_h2 + delta_p_f_h3
    delta_p_f_c = delta_p_f_c1 + delta_p_f_c2 + delta_p_f_c3

    #v_cooling = massRateCooling / (coolingDensity * (0.25 * math.pi * (port_D ** 2)))

    delta_p_p_h = 1.5 * (G_port_h ** 2) / (2 * PSI("D", "P", c20.p.val, "T", 273.15 + (c20.T.val + c21.T.val) / 2, coolingFluid2))
    delta_p_p_c = 1.5 * (G_port_c ** 2) / ((PSI("D", "P", pressureOut, "Q", 1, wf) + PSI("D", "P", pressureOut, "Q", 0, wf)))

    delta_p_c = delta_p_f_c + delta_p_p_c
    delta_p_h = delta_p_f_h + delta_p_p_h

    #Não se considera as perdas de carga por ação gravítica por causa da reduzida variação da mesma

    #Dimensionamento tubo alhetado

    air_velocity = 10 #m/s
    d_out = 0.01905
    tube_wall = 0.0015
    d_in = d_out - 2 * tube_wall
    r1_fin = d_out / 2
    L_fin = 0.00094  
    r2_fin = r1_fin + L_fin
    t_fin = 0.000279
    k_fin = 14.9
    k_tube = 14.9
    Lc_fin = L_fin + t_fin / 2
    Ap_fin = Lc_fin * t_fin
    A_fin = 2 * math.pi * ((r2_fin ** 2) - (r1_fin ** 2))
    r2c_fin = r2_fin + t_fin
    fins_per_inch = 28
    S_fin = 0.0254 / fins_per_inch - t_fin

    Re_outside = air_velocity * d_out * PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") / PSI("V", "T", 273.15 + T_amb, "P", 1.01325e5, "Air")

    if Re_outside >= 0.4 and Re_outside < 4:
        c = 0.989
        m_re = 0.33
    elif Re_outside >= 4 and Re_outside < 40:
        c = 0.911
        m_re = 0.385
    elif Re_outside >= 40 and Re_outside < 4000:
        c = 0.683
        m_re = 0.466
    elif Re_outside >= 4000 and Re_outside < 40000:
        c = 0.193
        m_re = 0.618
    else:
        c = 0.027
        m_re = 0.805

    Pr_outside = PSI("C", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * PSI("V", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") / PSI("L", "T", 273.15 + T_amb, "P", 1.01325e5, "Air")
    Nu_outside = c * (Re_outside ** m_re) * (Pr_outside ** (1/3))
    h_outside = Nu_outside * PSI("L", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") / d_out

    Pr_1 = PSI("D", "T", 273.15 + (c1.T.val + T_satWF) / 2, "P", pressureOut, wf) * PSI("V", "T", 273.15 + (c1.T.val + T_satWF) / 2, "P", pressureOut, wf) / PSI("L", "T", 273.15 + (c1.T.val + T_satWF) / 2, "P", pressureOut, wf)
    Re_stage1 = 4 * massRateWF / (math.pi * d_in * PSI("V", "T", 273.15 + (c1.T.val + T_satWF) / 2, "P", pressureOut, wf))
    Nu_stage1 = 0.023 * (Re_stage1 ** 0.8) * (Pr_1 ** 0.3)
    h_stage1 = Nu_stage1 * PSI("L", "T", 273.15 + (T_satWF + c1.T.val) / 2, "P", pressureOut, wf) / d_in

    Pr_3 = PSI("D", "T", 273.15 + (c2.T.val + T_satWF) / 2, "P", pressureOut, wf) * PSI("V", "T", 273.15 + (c2.T.val + T_satWF) / 2, "P", pressureOut, wf) / PSI("L", "T", 273.15 + (c2.T.val + T_satWF) / 2, "P", pressureOut, wf)
    Re_stage3 = 4 * massRateWF / (math.pi * d_in * PSI("V", "T", 273.15 + (c2.T.val + T_satWF) / 2, "P", pressureOut, wf))
    Nu_stage3 = 0.023 * (Re_stage3 ** 0.8) * (Pr_3 ** 0.3)
    h_stage3 = Nu_stage3 * PSI("L", "T", 273.15 + (T_satWF + c2.T.val) / 2, "P", pressureOut, wf) / d_in

    m = math.sqrt(2 * h_outside / (k_fin * t_fin))
    mr1 = m * r1_fin
    mr2 = m * r2c_fin
    I0_mr1 = iv(0, mr1)
    I0_mr2 = iv(0,mr2)
    I1_mr1 = iv(1, mr1)
    I1_mr2 = iv(1, mr2)
    K0_mr1 = kv(0, mr1)
    K0_mr2 = kv(0, mr2)
    K1_mr1 = kv(1, mr1)
    K1_mr2 = kv(1, mr2)

    eff_fin = (2 * r1_fin * (K1_mr1 * I1_mr2 - I1_mr1 * K1_mr2)) / (m *((r2c_fin ** 2) - (r1_fin ** 2)) * (K0_mr1 * I1_mr2 + I0_mr1 * K1_mr2))
    xx = 0.5
    Re_tube_l = 4 * massRateWF * xx / (math.pi * d_in * PSI("V", "T", 273.15 + T_satWF, "Q", 0, wf))
    Pr_tube_l = PSI("D", "T", 273.15 + T_satWF, "Q", 0, wf) * PSI("V", "T", 273.15 + T_satWF, "Q", 0, wf) / PSI("L", "T", 273.15 + T_satWF, "Q", 0, wf)
    Martinelli = (((1 - xx) / xx) ** 0.9) * ((PSI("D", "T", 273.15 + T_satWF, "Q", 1, wf) / PSI("D", "T", 273.15 + T_satWF, "Q", 0, wf)) ** 0.5) * ((PSI("V", "T", 273.15 + T_satWF, "Q", 0, wf) / PSI("V", "T", 273.15 + T_satWF, "Q", 1, wf)) ** 0.1)
    Nu_stage2 = 0.023 * (Re_tube_l ** 0.8) * (Pr_tube_l ** 0.4) * (1 + 2.22 / (Martinelli ** 0.89))
    h_stage2 = Nu_stage2 * PSI("L", "T", 273.15 + T_satWF, "Q", 0, wf) / d_in

    Q_stage1_L = (c1.T.val - T_amb) / ((1 / (2 * math.pi * r1_fin * h_stage1)) + (math.log(r2_fin / r1_fin) / (2 * math.pi * k_tube)) + ((S_fin + t_fin) / (h_outside * (2 * math.pi * r2_fin * S_fin + A_fin * eff_fin))))
    length_1 = ((c1.h.val * 1000 - PSI("H", "P", pressureOut, "Q", 1, wf)) * massRateWF) / Q_stage1_L

    Q_stage2_L = (T_satWF - T_amb) / ((1 / (2 * math.pi * r1_fin * h_stage2)) + (math.log(r2_fin / r1_fin) / (2 * math.pi * k_tube)) + ((S_fin + t_fin) / (h_outside * (2 * math.pi * r2_fin * S_fin + A_fin * eff_fin))))
    length_2 = (PSI("H", "P", pressureOut, "Q", 1, wf) - PSI("H", "P", pressureOut, "Q", 0, wf)) * massRateWF / Q_stage2_L

    Q_stage3_L = (T_satWF - T_amb) / ((1 / (2 * math.pi * r1_fin * h_stage3)) + (math.log(r2_fin / r1_fin) / (2 * math.pi * k_tube)) + ((S_fin + t_fin) / (h_outside * (2 * math.pi * r2_fin * S_fin + A_fin * eff_fin))))
    length_3 = ((PSI("H", "P", pressureOut, "Q", 0, wf) - c2.h.val * 1000) * massRateWF) / Q_stage3_L

    length_total = length_1 + length_2 + length_3

    Q_total_tube = massRateWF * (c1.h.val - c2.h.val) * 1000 #W
    Tf_air = T_amb + Q_total_tube / (massRateCooling1 * PSI("C", "T", 273.15 + T_amb, "P", 1.01325e5, "Air"))
    
    t_air_1 = Tf_air
    t_pack_1 = t_pack

    t0 = 0
    t1 = 1
    t_air_2 = t_pack_1 + (t_air_1 - t_pack_1) / (math.exp((transfer_area_pack * h_outside) / (massRateCooling1 * PSI("C", "T", 273.15 + t_air_1, "P", 1.01325e5, "Air"))))
    pack_heat_rate = massRateCooling1 * PSI("C", "T", 273.15 + t_air_1, "P", 1.01325e5, "Air") * (t_air_1 - t_air_2)
    t_pack_2 = t_pack_1 + pack_heat_rate * (t1 - t0) / (m_pack * cp_pack)
    t_pack_1 = t_pack_2

    while t_pack_2 < t_pack_final:
        t_air_2 = t_pack_1 + (t_air_1 - t_pack_1) / (math.exp((transfer_area_pack * h_outside) / (massRateCooling1 * PSI("C", "T", 273.15 + t_air_1, "P", 1.01325e5, "Air"))))
        pack_heat_rate = massRateCooling1 * PSI("C", "T", 273.15 + t_air_1, "P", 1.01325e5, "Air") * (t_air_1 - t_air_2)
        t_pack_2 = t_pack_1 + pack_heat_rate * (t1 - t0) / (m_pack * cp_pack)
        t_pack_1 = t_pack_2
        t0 += 1
        t1 += 1


    print("Área necessária (permutador de placas): " + str(round(area_needed, 3)) + " m2")
    print("Área catalogo (permutador de placas): " + str(round(area_catalogo,3)) + " m2")
    print("Perda de carga no lado da água: " + str(round(delta_p_c / 1000, 4)) + " kPa")
    print("Perda de carga no lado do fluído de trabalho: " + str(round(delta_p_h / 1000, 3)) + " kPa")
    print("O tubo alhetado tem um comprimento total de: " + str(round(length_total, 2)) + " m")
    print("O ar exterior sairá a uma temperatura média de " + str(round(Tf_air, 2)) + " ºC")
    print("Para que o pack de baterias atinja uma temperatura final de " + str(t_pack_final) + " ºC, sendo que começa a uma temperatura incial de " + str(round(t_pack, 1)) + " ºC, irá demorar aproximadamente " + str(round(t1 / 60)) + " minutos.")
    
    print("\n")

    print("Modo arrefecimento")

    heatPump_cooling = Network(T_unit='C', p_unit='Pa', h_unit='kJ / kg', iterinfo=False)

    cC0 = Connection(Cooling1, 'out2', compressor, 'in1', label= 'Estado inicial')
    cC1 = Connection(compressor, 'out1', Cooling2, 'in1', label= 'Pós compressão')
    cC2 = Connection(Cooling2, "out1", valve, "in1", label= "Pós arrefecimento") 
    cC3 = Connection(valve, "out1", closer, "in1", label= "Pós expansãp")
    cC4 = Connection(closer, "out1", Cooling1, "in2", label= "Cycle Closer")

    cC10 = Connection(source1, "out1", Cooling1, "in1", label= "Ambiente")
    cC11 = Connection(Cooling1, "out1", sink1, "in1", label= "Ambiente arrefecido")

    cC20 = Connection(source2, "out1", Cooling2, "in2", label= "Água fria")
    cC21 = Connection(Cooling2, "out2", sink2, "in1", label= "Água quente")

    heatPump_cooling.add_conns(cC0, cC1, cC2, cC3, cC4, cC10, cC11, cC20, cC21)

    wf = "R32"
    massRateWF = 0.037 #kg/s
    T_amb = 10
    T_cooling_in2 = 5 #Água
    coolingHeated1 = -5 #Ar
    coolingHeated2 = 8 #Água
    coolingVolumeRate2 = 2 / 3600 #m3/s
    coolingFluid2 = "Water"
    coolingFluid1 = "Air"
    coolingDensity2 = PSI("D", "T", 273.15 + T_cooling_in2, "P", 3e5, coolingFluid2) #kg/m3 Água
    massRateCooling2 = coolingVolumeRate2 * coolingDensity2 #kg/s
    coolingVolume1 = 800 / 3600 #m3/s
    massRateCooling1 = PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * coolingVolume1

    t_pack = T_amb
    t_pack_final = 5

    Cooling1.set_attr(pr1 = condenserPressureLossRatioHot, pr2 = 1)
    compressor.set_attr(eta_s = efficiencyCompressor, pr = pressureOut / pressureIn)
    Cooling2.set_attr(pr1 = evaporatorPressureLossRatioHot, pr2 = 1)

    cC0.set_attr(fluid={wf: 1}, m = massRateWF, T = T_wf_in, p = pressureIn)
    

    cC10.set_attr(fluid={coolingFluid1: 1}, p = 1.01325e5, T = T_amb, m = massRateCooling1)

    cC20.set_attr(fluid={coolingFluid2: 1}, p = 3e5, T = T_cooling_in2, m = massRateCooling2)


    valve.set_attr(pr = pressureIn / pressureOut)

    heatPump_cooling.solve(mode= 'design')
    heatPump_cooling.save('hp_design_cooling')
    heatPump_cooling.print_results()

Answered by fwitte

Aug 1, 2024

Hm... well it is kind of a bug, and kind of not. The heat exchanger cold side effectiveness uses the cold side pressure and the hot side inlet temperature of a heat exchanger (hot side effectiveness is hot side pressure and cold side inlet temperature). In this case, the temperature at the hot side inlet is out of bounds for the cold side fluid, which is unfortunate. I think it would be best, to except this issue in the postprocessing. I will do a quick fix for this this weekend. If you do not want to wait, I would like to invite you to make the change and open a pull request :)

Thank you for reporting and have a nice evening

View full answer

rafaelmarqferreira · 2024-07-19T08:52:43Z

rafaelmarqferreira
Jul 19, 2024
Author

PS: Trying to build a similar model without a cycle close seems to somewhat work, except for the fact that the end value won´t match the starting value

from tespy.networks import Network
from CoolProp.CoolProp import PropsSI as PSI
import matplotlib.pyplot as plt
from fluprodia import FluidPropertyDiagram
import math
from scipy import constants
import sys
from scipy.special import iv, kv
from tespy.connections import Connection
from tespy.components import (CycleCloser, Compressor, Valve, Pump, HeatExchanger, Source, Sink)

heatPump_cooling = Network(T_unit='C', p_unit='Pa', h_unit='kJ / kg', iterinfo=False)

closer = CycleCloser('cycle closer') #Chamar componente ciclo de fecho
Cooling1 = HeatExchanger('Cooling 1')
Cooling2 = HeatExchanger('Cooling 2') #Chamar componente permutador de calor
valve = Valve('expansion valve') #Chamar componente válvula
compressor = Compressor('compressor') #Chamar componente compressor
pump = Pump('pump') #Chamar componente bomba

source1 = Source('ambIn1')
source2 = Source('ambIn2')
sink1 = Sink('ambOut1')
sink2 = Sink('ambOut2')

test_source1 = Source('test_source1')
test_source2 = Source('test_source2')
test_source3 = Source('test_source3')
test_sink1 = Sink('test_sink1')
test_sink2 = Sink('test_sink2')
test_sink3 = Sink('test_sink3')

wf = "R32"
massRateWF = 0.0183 #kg/s
T_wf_in = -33
pressureIn = 1.8e5
pressureOut = 16e5
T_amb = 10
T_cooling_in2 = 5 #Água
coolingHeated1 = 1 #Ar
coolingHeated2 = 8 #Água
coolingVolumeRate2 = 2 / 3600 #m3/s
coolingFluid2 = "Water"
coolingFluid1 = "Air"
coolingDensity2 = PSI("D", "T", 273.15 + T_cooling_in2, "P", 3e5, coolingFluid2) #kg/m3 Água
massRateCooling2 = coolingVolumeRate2 * coolingDensity2 #kg/s
coolingVolume1 = 800 / 3600 #m3/s
massRateCooling1 = PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * coolingVolume1

condenserPressureLossRatioHot = 1
evaporatorPressureLossRatioHot = 1

efficiencyCompressor = 0.85

efficiencyPump = 0.65

Cooling1.set_attr(pr1 = condenserPressureLossRatioHot, pr2 = 1)
compressor.set_attr(eta_s = efficiencyCompressor, pr = pressureOut / pressureIn)
Cooling2.set_attr(pr1 = evaporatorPressureLossRatioHot, pr2 = 1)

valve.set_attr(pr = pressureIn / pressureOut)


c0 = Connection(test_source1, 'out1', compressor, 'in1', label= 'Estado inicial')
c1 = Connection(compressor, 'out1', Cooling2, 'in1', label= 'Pós compressão')
c2 = Connection(Cooling2, 'out1', valve, 'in1', label= 'Pós arrefecimento')
c3 = Connection(valve, 'out1', Cooling1, 'in2', label= 'Pós expansão')
c4 = Connection(Cooling1, 'out2', test_sink1, 'in1', label= 'Pós aquecimento')

c10 = Connection(test_source2, "out1", Cooling2, "in2", label= "Água fria")
c11 = Connection(Cooling2, "out2", test_sink2, "in1", label= "Água quente")

c20 = Connection(test_source3, 'out1', Cooling1, 'in1', label= 'Ar quente')
c21 = Connection(Cooling1, 'out1', test_sink3, 'in1', label= 'Ar frio')

c0.set_attr(fluid={wf: 1}, m = massRateWF, T = T_wf_in, p = pressureIn)

c10.set_attr(fluid={coolingFluid2: 1}, m = massRateCooling2, T = T_cooling_in2, p = 3e5)
c11.set_attr(T = coolingHeated2)

c20.set_attr(fluid={coolingFluid1: 1}, m = massRateCooling1, T = T_amb, p = 1.01325e5)
c21.set_attr(T = coolingHeated1)

heatPump_cooling.add_conns(c0, c1, c2, c3, c4, c10, c11, c20, c21)
heatPump_cooling.solve(mode='design')
heatPump_cooling.print_results()

result_dict = {}
result_dict.update({Cooling2.label: Cooling2.get_plotting_data()[1]})
result_dict.update({compressor.label: compressor.get_plotting_data()[1]})
result_dict.update({Cooling1.label: Cooling1.get_plotting_data()[2]})
result_dict.update({valve.label: valve.get_plotting_data()[1]})

diagram = FluidPropertyDiagram(wf)
diagram.set_unit_system(T='°C', p='Pa', h='kJ/kg')

for key, data in result_dict.items():
    result_dict[key]['datapoints'] = diagram.calc_individual_isoline(**data)

diagram.calc_isolines()

fig, ax = plt.subplots(1, figsize=(16, 10))
diagram.draw_isolines(fig, ax, 'logph', x_min= c3.h.val - 350, x_max= c1.h.val + 100, y_min= 1e4  , y_max= c1.p.val * 10)

for key in result_dict.keys():
    datapoints = result_dict[key]['datapoints']
    ax.plot(datapoints['h'], datapoints['p'], color='#ff0000')
    ax.scatter(datapoints['h'][0], datapoints['p'][0], color='#ff0000')

fig.savefig(str(wf) + '_logph_cooling.svg')

0 replies

fwitte · 2024-07-22T07:23:39Z

fwitte
Jul 22, 2024
Maintainer

Hi @rafaelmarqferreira and thanks for reaching out! I will respond to you later today :)

9 replies

rafaelmarqferreira Jul 29, 2024
Author

Seems to be working!
Thank you for your time!

fwitte Jul 30, 2024
Maintainer

One more side note to your earlier version of the cooling mode. The issue there seemed to be mostly related to the parametrization of the model: You had mass flow and temperature difference on both the heat source and the heat sink side specified. Something like this would work.

closer = CycleCloser('cycle closer') #Chamar componente ciclo de fecho
Cooling1 = HeatExchanger('Cooling 1')
Cooling2 = HeatExchanger('Cooling 2') #Chamar componente permutador de calor
valve = Valve('expansion valve') #Chamar componente válvula
compressor = Compressor('compressor') #Chamar componente compressor
pump = Pump('pump') #Chamar componente bomba

source1 = Source('ambIn1')
source2 = Source('ambIn2')
sink1 = Sink('ambOut1')
sink2 = Sink('ambOut2')

test_source1 = Source('test_source1')
test_source2 = Source('test_source2')
test_source3 = Source('test_source3')
test_sink1 = Sink('test_sink1')
test_sink2 = Sink('test_sink2')
test_sink3 = Sink('test_sink3')

wf = "R32"
massRateWF = 0.0183 #kg/s
T_wf_in = -33
pressureIn = 1.8e5
pressureOut = 16e5
T_amb = 10
T_cooling_in2 = 5 #Água
coolingHeated1 = 1 #Ar
coolingHeated2 = 8 #Água
coolingVolumeRate2 = 2 / 3600 #m3/s
coolingFluid2 = "Water"
coolingFluid1 = "Air"
coolingDensity2 = PSI("D", "T", 273.15 + T_cooling_in2, "P", 3e5, coolingFluid2) #kg/m3 Água
massRateCooling2 = coolingVolumeRate2 * coolingDensity2 #kg/s
coolingVolume1 = 800 / 3600 #m3/s
massRateCooling1 = PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * coolingVolume1

condenserPressureLossRatioHot = 1
evaporatorPressureLossRatioHot = 1

efficiencyCompressor = 0.85

efficiencyPump = 0.65

Cooling1.set_attr(pr1 = condenserPressureLossRatioHot, pr2 = 1)
compressor.set_attr(eta_s = efficiencyCompressor)
Cooling2.set_attr(pr1 = evaporatorPressureLossRatioHot, pr2 = 1)

valve.set_attr()

c0 = Connection(closer, 'out1', compressor, 'in1', label= 'Estado inicial')
c1 = Connection(compressor, 'out1', Cooling2, 'in1', label= 'Pós compressão')
c2 = Connection(Cooling2, 'out1', valve, 'in1', label= 'Pós arrefecimento')
c3 = Connection(valve, 'out1', Cooling1, 'in2', label= 'Pós expansão')
c4 = Connection(Cooling1, 'out2', closer, 'in1', label= 'Pós aquecimento')

c10 = Connection(test_source2, "out1", Cooling2, "in2", label= "Água fria")
c11 = Connection(Cooling2, "out2", test_sink2, "in1", label= "Água quente")

c20 = Connection(test_source3, 'out1', Cooling1, 'in1', label= 'Ar quente')
c21 = Connection(Cooling1, 'out1', test_sink3, 'in1', label= 'Ar frio')

c0.set_attr(fluid={wf: 1}, m = massRateWF, T = coolingHeated1 - 5, x=1)
c2.set_attr(T=coolingHeated2 + 5, x=0)

c10.set_attr(fluid={coolingFluid2: 1}, T = T_cooling_in2, p = 3e5)
c11.set_attr(T = coolingHeated2)

c20.set_attr(fluid={coolingFluid1: 1}, T = T_amb, p = 1.01325e5)
c21.set_attr(T = coolingHeated1)

nw.add_conns(c0, c1, c2, c3, c4, c10, c11, c20, c21)
nw.solve(mode='design')
nw.print_results()

If you have those mass flow and temperature difference values from measurements, instead of specifying them to the model, you can use them as check, if the model calculated the missing values to match the measurements correctly.

rafaelmarqferreira Aug 1, 2024
Author

Also noticed, which could be a bug but I'm oblivious to whether or not that's the case, but updating the current TESPy version led to the same error showing up, but with a different temperature value being printed on the terminal

Traceback (most recent call last):
  File "c:\Users\rferreira\Desktop\Coding_test\HP\NGS_hp.py", line 103, in <module>
    heatPump.solve(mode='design')
  File "C:\Users\rferreira\AppData\Local\miniforge3\envs\tespy-env\Lib\site-packages\tespy\networks\network.py", line 1995, in solve
    self.postprocessing()
  File "C:\Users\rferreira\AppData\Local\miniforge3\envs\tespy-env\Lib\site-packages\tespy\networks\network.py", line 2460, in postprocessing
    self.process_components()
  File "C:\Users\rferreira\AppData\Local\miniforge3\envs\tespy-env\Lib\site-packages\tespy\networks\network.py", line 2486, in process_components
    cp.calc_parameters()
  File "C:\Users\rferreira\AppData\Local\miniforge3\envs\tespy-env\Lib\site-packages\tespy\components\heat_exchangers\base.py", line 1067, in calc_parameters
    / self.calc_dh_max_cold()
      ^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\rferreira\AppData\Local\miniforge3\envs\tespy-env\Lib\site-packages\tespy\components\heat_exchangers\base.py", line 722, in calc_dh_max_cold
    h_at_T_in_hot = h_mix_pT(
                    ^^^^^^^^^
  File "C:\Users\rferreira\AppData\Local\miniforge3\envs\tespy-env\Lib\site-packages\tespy\tools\fluid_properties\functions.py", line 128, in h_mix_pT       
    return pure_fluid["wrapper"].h_pT(p, T)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\rferreira\AppData\Local\miniforge3\envs\tespy-env\Lib\site-packages\tespy\tools\fluid_properties\wrappers.py", line 204, in h_pT
    self.AS.update(CP.PT_INPUTS, p, T)
  File "CoolProp\\AbstractState.pyx", line 102, in CoolProp.CoolProp.AbstractState.update
  File "CoolProp\\AbstractState.pyx", line 104, in CoolProp.CoolProp.AbstractState.update
ValueError: Your temperature 390.646303 is not between 173.150000 and 373.150000.

 from tespy.networks import Network
from CoolProp.CoolProp import PropsSI as PSI
import matplotlib.pyplot as plt
from fluprodia import FluidPropertyDiagram
import math
from scipy import constants
from scipy.special import iv, kv
from tespy.connections import Connection
from tespy.components import (CycleCloser, Compressor, Valve, Pump, HeatExchanger, Source, Sink)

heatPump = Network(T_unit='C', p_unit='Pa', h_unit='kJ / kg', iterinfo=False)

closer = CycleCloser('cycle closer') #Chamar componente ciclo de fecho
Cooling1 = HeatExchanger('Permutador de placas Heating / Tubo alhetado Cooling')
Cooling2 = HeatExchanger('Tubo alhetado Heating / Permutador de placas Cooling') #Chamar componente permutador de calor
valve = Valve('expansion valve') #Chamar componente válvula
compressor = Compressor('compressor') #Chamar componente compressor
pump = Pump('pump') #Chamar componente bomba

source1 = Source('ambIn1')
source2 = Source('ambIn2')
sink1 = Sink('ambOut1')
sink2 = Sink('ambOut2')


c0 = Connection(Cooling2, 'out2', compressor, 'in1', label= 'Estado inicial')
c1 = Connection(compressor, 'out1', Cooling1, 'in1', label= 'Pós compressão')
c2 = Connection(Cooling1, "out1", closer, "in1", label= "Pós arrefecimento")
c3 = Connection(closer, "out1", valve, "in1", label= "Cycle Closer")
c4 = Connection(valve, "out1", Cooling2, "in2", label= "Pós expansão")

c10 = Connection(source1, "out1", Cooling1, "in2", label= "Água")
c11 = Connection(Cooling1, "out2", sink1, "in1", label= "Água aquecida")

c20 = Connection(source2, "out1", Cooling2, "in1", label= "Ambiente antes")
c21 = Connection(Cooling2, "out1", sink2, "in1", label= "Ambiente depois")


heatPump.add_conns(c0, c1, c2, c3, c4, c10, c11, c20, c21)

wf = "R32"
T_wf_in = -34
pressureIn = 1.8e5
pressureOut = 17e5
h0 = PSI("H", "P", pressureIn, "T", 273.15 + T_wf_in, wf) / 1000
fluidCheck = PSI("T", "P", pressureIn, "Q", 1, wf) - 273.15
T_amb = -30
T_cooling_in1 = 6
T_cooling_in2 = T_amb
coolingHeated1 = 12
mix = 'INCOMP::MPG[0.6]'
coolingFluid1 = mix
# coolingFluid1 = "Water"
coolingFluid2 = "Air"
coolingVolumeRate1 = 2.3 / 3600 #m3/s
coolingVolumeRate2 = 2000 / 3600 #m3/s
coolingDensity1 = PSI("D", "T", 273.15 + T_cooling_in1, "P", 3e5, coolingFluid1) #kg/m3
massRateCooling1 = coolingVolumeRate1 * coolingDensity1 #kg/s
massRateCooling2 = PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * coolingVolumeRate2
# print(PSI("V", "P", 3e5, "T", T_cooling_in1 + 273.15, mix))

t_pack = T_amb
t_pack_final = 20
m_pack = 335
cp_pack = 500
a_pack = 1.073 #m
b_pack = 0.578 #m
c_pack = 0.586 #m
volume_pack = a_pack * b_pack * c_pack #m3
transfer_area_pack = 4 * a_pack * c_pack + a_pack * b_pack #m2

if T_wf_in <= fluidCheck:
    print("Existência de líquido / mistura bifásica no compressor. Aumente a temperatura de entrada.")
else:
    print("Modo de aquecimento")
    condenserPressureLossRatioHot = 1
    evaporatorPressureLossRatioHot = 1

    efficiencyCompressor = 0.85

    efficiencyPump = 0.65

    Cooling1.set_attr(pr1 = 0.992, pr2 = 0.901)
    compressor.set_attr(eta_s = efficiencyCompressor)
    Cooling2.set_attr(pr1 = evaporatorPressureLossRatioHot, pr2 = 1)

    valve.set_attr(pr = pressureIn / pressureOut)

    c0.set_attr(fluid={wf: 1}, T = T_wf_in, p = pressureIn, m = 0.018)
    c2.set_attr(T = 22)
    
    c10.set_attr(fluid={coolingFluid1: 1}, p = 3e5, T = T_cooling_in1)
    c11.set_attr(T = coolingHeated1)

    c20.set_attr(fluid={coolingFluid2: 1}, p = 1.01325e5, T = T_cooling_in2)
    c21.set_attr(m = massRateCooling2)
    
    h =  PSI("H", "T", 273.15 + T_cooling_in1, "P", 3e5, coolingFluid1) / 1e3
    [c.set_attr(h0=h0) for c in [c10]]
    h =  PSI("H", "T", 273.15 + coolingHeated1, "P", 3e5, coolingFluid1) / 1e3
    [c.set_attr(h0=h0) for c in [c11]]

    heatPump.solve(mode='design')

    heatPump.save('hp_design_heating')
    heatPump.print_results()

fwitte Aug 1, 2024
Maintainer

Hm... well it is kind of a bug, and kind of not. The heat exchanger cold side effectiveness uses the cold side pressure and the hot side inlet temperature of a heat exchanger (hot side effectiveness is hot side pressure and cold side inlet temperature). In this case, the temperature at the hot side inlet is out of bounds for the cold side fluid, which is unfortunate. I think it would be best, to except this issue in the postprocessing. I will do a quick fix for this this weekend. If you do not want to wait, I would like to invite you to make the change and open a pull request :)

Thank you for reporting and have a nice evening

Answer selected by fwitte

fwitte Aug 4, 2024
Maintainer

I released a fixed version!

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Heat pump issue (Cooling mode) #528

{{title}}

Replies: 2 comments 9 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

Select a reply

Heat pump issue (Cooling mode) #528

rafaelmarqferreira Jul 19, 2024

Replies: 2 comments · 9 replies

rafaelmarqferreira Jul 19, 2024 Author

fwitte Jul 22, 2024 Maintainer

rafaelmarqferreira Jul 29, 2024 Author

fwitte Jul 30, 2024 Maintainer

rafaelmarqferreira Aug 1, 2024 Author

fwitte Aug 1, 2024 Maintainer

fwitte Aug 4, 2024 Maintainer

rafaelmarqferreira
Jul 19, 2024

Replies: 2 comments 9 replies

rafaelmarqferreira
Jul 19, 2024
Author

fwitte
Jul 22, 2024
Maintainer

rafaelmarqferreira Jul 29, 2024
Author

fwitte Jul 30, 2024
Maintainer

rafaelmarqferreira Aug 1, 2024
Author

fwitte Aug 1, 2024
Maintainer

fwitte Aug 4, 2024
Maintainer