Sizing pressure losses with Kv values

[rafaelmarqferreira]

Still the same model that I've been working on. Trying to introduce pressure losses on valves while knowing the volumetric rate and as well the specific Kv for said valve. Was able to run a first simulation by setting the pressure ratio to 1 (as a starting point) to get the correct density at each connection and extract the correct volumetric rate. After setting the characteristic pressure loss curve, it seems as though there is a missing parameter. What could it be? What recommendations would you give to solve this case?
Also, in the same subject, if the pressure loss was to be a set value of X Pa (or any other unit), how should it be declared?

Best wishes and thankful for the help up until now

from tespy.networks import Network
from tespy.tools.characteristics import CharLine
from CoolProp.CoolProp import PropsSI as PSI
import matplotlib.pyplot as plt
from fluprodia import FluidPropertyDiagram
import math
import numpy as np
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') 
Cooling1 = HeatExchanger('Permutador de placas Heating / Tubo alhetado Cooling')
Cooling2 = HeatExchanger('Tubo alhetado Heating / Permutador de placas Cooling')
safety = Valve('Válvula segurança')
needle1 = Valve('Válvula agulha 1')
needle2 = Valve('Válvula agulha 2')
valve1 = Valve('expansion valve 1') 
valve2 = Valve('expansion valve 2') 
valve_three_wayA = Valve('Válvula 3 vias A')
valve_three_wayB = Valve('Válvula 3 vias B')
circuit_filter = Valve('Filtro')
compressor = Compressor('compressor') 
pump = Pump('pump') 

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

c0 = Connection(circuit_filter, 'out1', compressor, 'in1', label= 'Estado inicial')
c1 = Connection(compressor, 'out1', safety, 'in1', label= 'Pós compressão')
c02 = Connection(safety, 'out1', valve_three_wayA, 'in1', label= 'Válvula segurança')
c03 = Connection(valve_three_wayA, 'out1', needle1, 'in1', label= 'Saída válvula 3 vias')
c04 = Connection(needle1, 'out1', Cooling1, 'in1', label= 'Válvula agulha 1')
c2 = Connection(Cooling1, "out1", closer, "in1", label= "Pós arrefecimento")
c3 = Connection(closer, "out1", needle2, "in1", label= "Cycle Closer")
c05 = Connection(needle2, 'out1', valve1, 'in1', label= 'Válvula agulha 2')
c31 = Connection(valve1, "out1", valve2, "in1", label= "Pós expansão inicial")
c4 = Connection(valve2, "out1", Cooling2, "in2", label= "Pós expansão final")
c00 = Connection(Cooling2, 'out2', valve_three_wayB, 'in1', label= 'Retorno válvula 3 vias')
c01 = Connection(valve_three_wayB, 'out1', circuit_filter, 'in1', label= 'Passagem pelo filtro')

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, c02, c03, c04, c2, c3, c05, c31, c4, c00, c01, c10, c11, c20, c21)

wf = "R32"
massRateWF = 0.035
T_wf_in = -34
pressureIn = 2e5
pressureOut = 12.5e5
# volumeRateWF = massRateWF / PSI("D", "T", 273.15 + c1.T.val, "P", pressureOut, wf) * 3600
pressureMix = 5e5
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", pressureMix, coolingFluid1) #kg/m3
massRateCooling1 = coolingVolumeRate1 * coolingDensity1 #kg/s
massRateCooling2 = PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * coolingVolumeRate2

# mass_flow_needle = np.array([0, 0.5, 1, 1.5, 2])
# pressure_drop_needle = np.array([0, 36762.5, 73525, 110287.5, 147050])
# kv_char = CharLine(x=mass_flow_needle, y=pressure_drop_needle)



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)

    c0.set_attr(fluid={wf: 1}, T = T_wf_in, p = pressureIn, m = massRateWF)
    c2.set_attr(x = 0)
    c31.set_attr(x = 0.1)
    safety.set_attr(pr = 1)
    valve_three_wayA.set_attr(pr = 1)
    valve_three_wayB.set_attr(pr = 1)    
    needle1.set_attr(pr = 1)
    needle2.set_attr(pr = 1)
    circuit_filter.set_attr(pr = 1)

    c10.set_attr(fluid={coolingFluid1: 1}, p = pressureMix, 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", pressureMix, coolingFluid1) / 1e3
    [c.set_attr(h0=h0) for c in [c10, c11]]

    heatPump.solve(mode='design')

    density = PSI("D", "P", c1.p.val, "T", 273.15 + c1.T.val, wf)

    volu_flow_safety = np.array([0 / density, 0.5 / density, 1 / density, 1.5 / density, 2 / density])
    pressure_drop_safety = np.array([0, 2.75e5, 5.5e5, 8.25e5, 11e5])
    kv_char_s = CharLine(x=volu_flow_safety, y=pressure_drop_safety)

    safety.set_attr(pr = None, dp_char={'char_func': kv_char_s})

    heatPump.solve(mode='design')

[fwitte]

Hi!

You if you provide it like that, you have to add is_set

    safety.set_attr(pr = None, dp_char={'char_func': kv_char_s, "is_set": True})

The API could probably be designed in a better way (or better documented?) here, right now I do not have the capacities to tackle that. I'd invite you to do the changes. If you are interested, I could give you a short introduction on the set_attr method (which is quite a mess right now), and you could make a suggestion on how to improve it :).

Best

[fwitte]

By the way, you can easily integrate the pressure loss equation into the Valve (or any) component. So you would not need to precalculate pressure loss in advance and then insert the lookup table. Can you post the mathematical relationship here, I can show you how to integrate it.

[rafaelmarqferreira]

First of all, thank you for your continuous availability.
A work around that I've tried is to firstly assume that there's no pressure losses in the valve instances, so that proper starting values could be achieve. After that, I introduced the pressure losses on each valve 1 by 1 depending on the Kv values given on the components data sheet, as such in the following code. I overwrote the pr constant where it previously were set to 1.

from tespy.networks import Network
from tespy.tools.characteristics import CharLine
from CoolProp.CoolProp import PropsSI as PSI
import matplotlib.pyplot as plt
from fluprodia import FluidPropertyDiagram
import math
import numpy as np
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') 
Cooling1 = HeatExchanger('Permutador de placas Heating / Tubo alhetado Cooling')
Cooling2 = HeatExchanger('Tubo alhetado Heating / Permutador de placas Cooling')
safety = Valve('Válvula segurança')
needle1 = Valve('Válvula agulha 1')
needle2 = Valve('Válvula agulha 2')
valve1 = Valve('expansion valve 1') 
valve2 = Valve('expansion valve 2') 
valve_three_wayA = Valve('Válvula 3 vias A')
valve_three_wayB = Valve('Válvula 3 vias B')
circuit_filter = Valve('Filtro')
compressor = Compressor('compressor') 
pump = Pump('pump') 

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

c0 = Connection(circuit_filter, 'out1', compressor, 'in1', label= 'Estado inicial')
c1 = Connection(compressor, 'out1', safety, 'in1', label= 'Pós compressão')
c02 = Connection(safety, 'out1', valve_three_wayA, 'in1', label= 'Válvula segurança')
c03 = Connection(valve_three_wayA, 'out1', needle1, 'in1', label= 'Saída válvula 3 vias')
c04 = Connection(needle1, 'out1', Cooling1, 'in1', label= 'Válvula agulha 1')
c2 = Connection(Cooling1, "out1", closer, "in1", label= "Pós arrefecimento")
c3 = Connection(closer, "out1", needle2, "in1", label= "Cycle Closer")
c05 = Connection(needle2, 'out1', valve1, 'in1', label= 'Válvula agulha 2')
c31 = Connection(valve1, "out1", valve2, "in1", label= "Pós expansão inicial")
c4 = Connection(valve2, "out1", Cooling2, "in2", label= "Pós expansão final")
c00 = Connection(Cooling2, 'out2', valve_three_wayB, 'in1', label= 'Saída permutador ar')
c01 = Connection(valve_three_wayB, 'out1', circuit_filter, 'in1', label= 'Saída válvula 3 vias / Entrada no filtro')

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, c02, c03, c04, c2, c3, c05, c31, c4, c00, c01, c10, c11, c20, c21)

wf = "R32"
massRateWF = 0.035
T_wf_in = -34
pressureIn = 1.8e5
pressureOut = 12.5e5
# volumeRateWF = massRateWF / PSI("D", "T", 273.15 + c1.T.val, "P", pressureOut, wf) * 3600
pressureMix = 5e5
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", pressureMix, coolingFluid1) #kg/m3
massRateCooling1 = coolingVolumeRate1 * coolingDensity1 #kg/s
massRateCooling2 = PSI("D", "T", 273.15 + T_amb, "P", 1.01325e5, "Air") * coolingVolumeRate2

kv_safety = 5.5
kv_needle1 = 0.85 * 0.865
kv_needle2 = 0.85 * 0.865
delta_p_filter = 0.21e5
delta_p_3way = 0.1e5

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)

    c0.set_attr(fluid={wf: 1}, T = T_wf_in, p = pressureIn, m = massRateWF)
    c2.set_attr(x = 0)
    c31.set_attr(x = 0.01)
    safety.set_attr(pr = 1)
    valve_three_wayA.set_attr(pr = 1)
    valve_three_wayB.set_attr(pr = 1)    
    needle1.set_attr(pr = 1)
    needle2.set_attr(pr = 1)
    circuit_filter.set_attr(pr = 1)

    c10.set_attr(fluid={coolingFluid1: 1}, p = pressureMix, 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", pressureMix, coolingFluid1) / 1e3
    [c.set_attr(h0=h0) for c in [c10, c11]]

    heatPump.solve(mode='design')

    density_safety = PSI("D", "P", c1.p.val, "T", 273.15 + c1.T.val, wf)
    safety.set_attr(pr = 1 - (kv_safety * (3600 * massRateWF / density_safety)) / pressureOut)

    heatPump.solve(mode='design')

    valve_three_wayA.set_attr(pr = (c02.p.val - delta_p_3way) / c02.p.val)

    heatPump.solve(mode='design') 

    density_needle1 = PSI("D", "P", c03.p.val, "T", 273.15 + c03.T.val, wf)
    needle1.set_attr(pr = 1 - (kv_needle1 * (3600 * massRateWF / density_needle1)) / c03.p.val)

    heatPump.solve(mode='design') 
    
    density_needle2 = PSI("D", "P", c2.p.val, "Q", c2.x.val, wf)
    needle2.set_attr(pr = 1 - (kv_needle2 * (3600 * massRateWF / density_needle2)) / c2.p.val)

    heatPump.solve(mode='design') 

    valve_three_wayB.set_attr(pr = (c00.p.val - delta_p_3way) / c00.p.val)

    heatPump.solve(mode='design') 

    circuit_filter.set_attr(pr = (c01.p.val - delta_p_filter) / c01.p.val)

    heatPump.solve(mode='design')

Given the Kv coefficient [m3/h] of a valve, the pressure loss in bar would be given as delta_p = ((volumetric_flow/Kv) ** 2). Having said that, the equation on the bottom can be derived to rewrite the pr value. Whenever a flat pressure was assumed (lack of Kv for now), the top equation can be used (apologies for the mess on the picture).

As for valve1 and valve2, those represent the same valve in essence, since I've tried to integrate the length needed for a capillary tube so that the desired outlet pressure was achieved. I'll leave the part of the script here:

d_cap = 6.35e-3
    a_cap = math.pi * (d_cap ** 2) / 4
    delta_L = 0
    t1 = c2.T.val
    p1 = c2.p.val
    v1 = 4 * massRateWF / (math.pi * PSI("D", "P", c2.p.val, "H", 1000 * c2.h.val, wf) * (d_cap ** 2))
    t2 = t1 - 0.01
    c31.set_attr(x = None, T = t2)
    heatPump.solve(mode='design')
    x = c31.x.val
    p2 = c31.p.val
    v2 =  4 * massRateWF / (math.pi * PSI("D", "T", 273.15 + t2, "Q", c31.x.val, wf) * (d_cap ** 2))
    vm = (v1 + v2) / 2
    delta_L_inc = ((p1 - p2) * (a_cap / massRateWF) + (v1 - v2)) * (2 * d_cap / vm)
    delta_L = delta_L + delta_L_inc
    while abs(p2 - c4.p.val) / c4.p.val >= 0.01:
        t1 = t2
        t2 = t1 - 0.1
        p1 = p2
        v1 = v2
        c31.set_attr(x = None, T = t2)
        heatPump.solve(mode='design')
        x = c31.x.val
        p2 = c31.p.val
        v2 =  4 * massRateWF / (math.pi * PSI("D", "T", 273.15 + t2, "Q", c31.x.val, wf) * (d_cap ** 2))
        vm = (v1 + v2) / 2
        delta_L_inc = ((p1 - p2) * (a_cap / massRateWF) + (v1 - v2)) * (2 * d_cap / vm)
        delta_L = delta_L + delta_L_inc
        
        

Hopefully I was clear enough! Thank you for the continuous support!
Best