Compute an exact pressure loss in a vessel

[rafaelmarqferreira]

So i have the following P-H diagram for R32 whose connection c2 (top-left point) has a pressure value of 13 bar and a temperature value of 8 ºC. Is there anyway to compute the specific pressure value at which, in a capillary tube, the refrigerant will be at a quality value of 0, at the same enthalpy value?

Calling the CoolProp library as such

PSI("P", "H", 1000 * c2.h.val, "Q", 0, wf)

will lead to the following error: ValueError: Input pair variable is invalid and output(s) are non-trivial; cannot do state update : PropsSI("P","H",214140.6804,"Q",0,"R32")

The goal for this would be to model the necessary length for the capillary tube to achieve a pressure loss of 11,2 bar (outlet of 1,8 bar), starting from a subcooled fluid to a two-stage refrigerant

[rafaelmarqferreira]

After testing, tried using 2 back to back expansion valves, which apparently does the job fairly well. Guess I'll just leave the code as an example for futures cases. Sorry for the unnecessary thread but oh well, at least it might help someone else in the future

from tespy.networks import Network
from CoolProp.CoolProp import PropsSI as PSI
import matplotlib.pyplot as plt
from fluprodia import FluidPropertyDiagram
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')
valve1 = Valve('expansion valve 1') #Chamar componente válvula
valve2 = Valve('expansion valve 2') #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", valve1, "in1", label= "Cycle Closer")
c4 = Connection(valve1, "out1", valve2, "in1", label= "Pós expansão inicial")
c5 = Connection(valve2, "out1", Cooling2, "in2", label= "Pós expansão final")

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, c5, c10, c11, c20, c21)

wf = "R32"
massRateWF = 0.01
T_wf_in = -34
pressureIn = 1.8e5
pressureOut = 13e5
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 = 0
T_cooling_in2 = T_amb
coolingHeated1 = 8
mix = 'INCOMP::GKN[0.6]'
coolingFluid1 = mix
coolingFluid2 = "Air"
coolingVolumeRate1 = 2.8 / 3600 #m3/s
coolingVolumeRate2 = 20000 / 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

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 = 1, pr2 = 1)
    compressor.set_attr(eta_s = efficiencyCompressor, pr = pressureOut / pressureIn)
    Cooling2.set_attr(pr1 = 1, pr2 = 1)

    # valve.set_attr(pr = pressureIn / pressureOut)

    c0.set_attr(fluid={wf: 1}, T = T_wf_in, p = pressureIn, m = massRateWF)
    c2.set_attr(T = 8)
    c4.set_attr(x = 0)
    # c5.set_attr(p = pressureIn)
    
    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, c11]]

    heatPump.solve(mode='design')
    # heatPump.save('hp_design_heating')
    # heatPump.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({valve1.label: valve1.get_plotting_data()[1]})
    result_dict.update({valve2.label: valve2.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_heating.svg')
    heatPump.print_results()

[fwitte]

Hi @rafaelmarqferreira and thank you for updating.

I think you could also do a simple Newton search instead:

from CoolProp.CoolProp import PropsSI


h_target = 214140.6804
p = 13e5
d = 1e-3

while True:
    residual = h_target - PropsSI("H", "P", p, "Q", 0, "R32")
    derivative = ((-PropsSI("H", "P", p + d, "Q", 0, "R32")) - (-PropsSI("H", "P", p - d, "Q", 0, "R32"))) / (2 * d)
    p -= residual / derivative
    if abs(residual) < 1e-3:
        break

print(p)

1042459.6668

Best!