In [40]:

    

#from __future__ import division, print_function
import math
import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
import CoolProp.CoolProp as CP

## Pipe dimensions
L = 10.0
dy_mm = 19.05 
z_wall_mm = 0.81
di_mm = dy_mm - 2*z_wall_mm 
z_wall = z_wall_mm / 1000
di =di_mm/1000  
e_mm = 0.0015  #for new copper pipes e is [0.0013,0.0015]
e= e_mm/1000
z_iso =  10/1000 #thickness-insulation
λambda_iso = 0.039 #conductivity-insulation

#Inputs 
Q_dot_cool = np.arange(5e3, 15e3, 100)
#Q_dot_cool = 10e3
dT_sh = 5
T_dew = 273.15 - 10
T_i = T_dew + dT_sh 
T_A = 273.15 + 25
T_c = 30 + 273.15
dT_sc = 2

old_fluids = ['R404A','R290']
fluids = [
    {'name': 'R404A', 'color': 'c--'},
    {'name': 'R290', 'color': 'r--'},
    {'name': 'R600a', 'color': 'g--'}
]

################Gas Pipes####################

ax1 = plt.subplot(2,1,1)
ax1.set_xlabel('Cooling capacity / kW')
ax1.set_ylabel('Pressure drop / kPa')

ax2 = plt.subplot(212, sharex=ax1)
ax2.set_xlabel('Cooling capacity / kW')
ax2.set_ylabel('Heat ingress / W')

for dict_fluid in fluids:
    
    fluid = dict_fluid['name']
    color = dict_fluid['color']
    
    p_c = CP.PropsSI('P','Q', 0, 'T', T_c, fluid)
    p_i =CP.PropsSI('P','Q', 1, 'T', T_dew, fluid)

    ## Fluid properties
    rho_i = CP.PropsSI('D','T', T_i, 'P', p_i, fluid)
    cp = CP.PropsSI('C','P', p_i, 'T', T_i, fluid) 
    mu = CP.PropsSI('V','T', T_i, 'P', p_i, fluid) #dynamic viscocity
    nu = mu / rho_i #kinematic
    λambda = CP.PropsSI('L','T', T_i, 'P', p_i, fluid)
    Pr = cp*mu/λambda

    #dh_sat = CP.PropsSI("H", "Q", 1, 'P', p_i, fluid)
    h_out = CP.PropsSI('H','P', p_i, 'T', T_i, fluid)
    h_in = CP.PropsSI('H','P', p_c, 'T', T_c-dT_sc, fluid)
    dh_s = h_out - h_in

    #print("hin={0}, hout={1}, dh={2} ".format(h_in,h_out,dh_s))

    ## Flow calculations
    m_dot = Q_dot_cool / dh_s
    V_dot = m_dot / rho_i
    c_out = V_dot / ( math.pi * di*di /4.0 )
    Re_D = c_out*di/nu #kinematic visc

    #### Heat Transfer part ####

    if np.all(Re_D <= 2000):
        Re_D = 2000

    if np.all(Re_D > 0): 
            Nus_D = 0.0235*(Re_D**0.8-230)*(1.8*Pr**0.3-0.8) 
    else: 
            Nus_D = 1

    alpha_i = Nus_D*λambda/di
    alpha_o = 10 

    dy = di+2*z_wall
    diso = dy+2*z_iso
    A_i = math.pi*di*L
    R_i = 1/(A_i*alpha_i)

    R_iso = np.log(diso/dy)/(2*math.pi*L*λambda_iso)
    A_o = math.pi*diso*L
    R_o = 1/(A_o*alpha_o)

    R_TOT = R_i+R_iso+R_o

    Q = (T_A-T_i)/R_TOT
    T_int = Q/(m_dot*cp)+T_i

    T_surf = (T_int+T_i)/2+Q*(R_i+R_iso)

    ## Pressure Drop

    if np.all((Re_D<49000)): 
        f = 0.316*Re_D**(-0.25)
    if np.all((Re_D>=49000)): 
        f = 0.184*Re_D**(-0.2) 
    if np.all((Re_D<=2000)): 
        f = 64/Re_D
    else:
        eD=e/di
        f_0 = 0.25*(np.log10(eD/3.7+5.74*Re_D**(-0.9)))**(-2)
        f = (1/(-2*np.log10(eD/3.7+2.51/(Re_D*f_0**0.5))))**2

    while True:
        f_0=f   
        f = (1/(-2*np.log10(eD/3.7+2.51/(Re_D*f_0**0.5))))**2

        if np.all((abs(f_0-f)<0.000001)):
            break

    L_E=0
    DpDx = 0.5*f*rho_i*c_out*c_out/di
    DpF = DpDx*(L+L_E)
    DpE = L_E*0.5*rho_i*c_out*c_out
    Dp_g = DpF+DpE
    p_o = p_i-Dp_g
    
    ax1.plot(Q_dot_cool/1e3,Dp_g/1e3,color,label=fluid)
    ax2.plot(Q_dot_cool/1e3,Q,color,label=fluid)
    #plt.setp(ax2.get_xticklabels(), visible=False)
    #setp(ax2, xticklabels=[])
    

#print("di={0}, velo={1}, V_dot={2}, m_dot={3},Q={4},Re={5}, dens={6}, DPg={7}".format(di,c_out,V_dot*3600.0,m_dot,Q,Re_D,rho_i,Dp_g))

#line0,=plt.plot(Q_dot_cool,Dp_g,'c--',label="R404A")
#line1,=plt.plot(Q_dot_cool,Dp_g,'r--',label="R290")
#line2,=plt.plot(Q,Dp_g,'b--',label="R404A")

### CoolPack for pressure drops ###
ax1.plot(np.linspace(5, 15, 3),np.array([9.5,34.3,73.4]),'gx',label="R600a")
ax1.plot(np.linspace(5, 15, 3),np.array([6.2,22.6,48.7]),'cx',label="R404A")
ax1.plot(np.linspace(5, 15, 3),np.array([3.3,12,25.5]),'rx',label="R290")

### CoolPack for heat ingress ###
ax2.plot(np.linspace(5, 15, 3),np.array([78.71,79.31,79.53]),'gx',label="R600a")
ax2.plot(np.linspace(5, 15, 3),np.array([78.9,79.42,79.62]),'cx',label="R404A")
ax2.plot(np.linspace(5, 15, 3),np.array([78.7,79.3,79.5]),'rx',label="R290")


#plt.legend(bbox_to_anchor=(1, 1), loc=0, borderaxespad=0.5)
plt.sca(ax1)
plt.legend(bbox_to_anchor=(0., 1.05, 1., .105), loc=3,
           ncol=3, mode="expand", borderaxespad=0.)
#ax = plt.gca()

plt.tight_layout()


    

In [ ]:

    




    

In [ ]: