#!/usr/bin/env python
# coding: utf-8

# In[1]:


import math

import CoolProp
import CoolProp.CoolProp as CP


# ## Algorytm obliczeń współczynnika wnikania ciepła po stronie wrzącego czynnika
# Podczas wrzenia czynnika w rurkach poziomych.

# ### Założenia projektowe

# In[2]:


Q_k = 15 * 1e3 # [kW] - wydajność grzewcza projektowanej pompy ciepła

T_0 = 3 + 273.15 # [K] - temperatura wrzenia czynnika chłodniczego
T_k = 40 + 273.15 # [K] - temperatura skraplania czynnika chłodniczego
dT_sup = 10 # [°C] - przegrzanie

n_iz = 0.5 # założony współczynnik sprawności izentropowej sprężarki


# In[3]:


# czynnik chłodniczy
refrigerant = "R134a"


# ### Obliczanie punktów charakterystycznych obiegu pompy ciepła

# In[4]:


p_0 = CP.PropsSI('P','Q',1,'T',T_0,refrigerant)
p_k = CP.PropsSI('P','Q',1,'T',T_k,refrigerant)

if dT_sup > 0: h1 = CP.PropsSI('H','P',p_0,'T',T_0+dT_sup,refrigerant)
else: h1 = CP.PropsSI('H','Q',1,'T',T_0,refrigerant)
s1 = CP.PropsSI('S','H',h1,'P',p_0,refrigerant)
h2 = CP.PropsSI('H','S',s1,'P',p_k,refrigerant)
h3 = CP.PropsSI('H','Q',0,'T',T_k,refrigerant)
h4 = h3

T_2 = CP.PropsSI('T','S',s1,'P',p_k,refrigerant)

# Korekta położenia punktu 2 po uwzględnieniu sprawności izentropowej sprężarki
h2b = h1 + (h2 - h1) / n_iz

q_0 = h1 - h4
q_k = h2b - h3

# strumień masy czynnika chłodniczego
m_0 = Q_k / q_k
print ("* Strumień masy czynnika chłodniczego: %.3f [kg/s]" % m_0)

Q_0 = m_0 * q_0
print ("* Wydajność chłodnicza parownika: %.3f [kW]" % (Q_0/1e3))


# ### Dane podstawowe

# In[5]:


T_p1 = 12 + 273.15 # [K] - temperatura powietrza na wlocie do parownika
T_p2 = 7 + 273.15 # [K] - temperatura powietrza na wylocie z parownika

# Średnia logarytmiczna różnica temperatur wymiennika
LMTD = (T_p1 - T_0) - (T_p2 - T_0) / math.log((T_p1 - T_0)/(T_p2 - T_0))

# Współczynnik wnikania ciepła po stronie powietrza (powinien być wcześniej obliczony)
α_p = 70 # [W/m2K]


# ### Czynnik chłodniczy i jego parametry w założonych temperaturach pracy wymiennika

# In[6]:


# ciśnienie nasycenia czynnika chłodniczego
p_n0 = CP.PropsSI('P','Q',0,'T',T_0,refrigerant)
print ("* Ciśnienie nasycenia %s: %.2f [J/kgK]" % (refrigerant,p_n0))

# ciepło właściwe czynnika chłodniczego
c_p0l = CP.PropsSI('C','Q',0,'T',T_0,refrigerant)
print ("* Ciepło właściwe %s (ciecz): %.2f [J/kgK]" % (refrigerant,c_p0l))
c_p0v = CP.PropsSI('C','Q',1,'T',T_0,refrigerant)
print ("* Ciepło właściwe %s (gaz): %.2f [J/kgK]" % (refrigerant,c_p0v))


# ciepło przemiany fazowej
r_0 = CP.PropsSI('H','Q',1.0,'T',T_0,refrigerant) - CP.PropsSI('H','Q',0,'T',T_0,refrigerant)
print ("* Ciepło przemiany fazowej %s: %.2f [J/kg]" % (refrigerant,r_0))

# gęstość czynnika chłodniczego (ciecz i gaz)
ρ_0l = CP.PropsSI('D','Q',0,'T',T_0,refrigerant)
print ("* Gęstość %s (ciecz): %.2f [kg/m3]" % (refrigerant,ρ_0l))
ρ_0v = CP.PropsSI('D','Q',1,'T',T_0,refrigerant)
print ("* Gęstość %s (gaz): %.2f [kg/m3]" % (refrigerant,ρ_0v))

# lepkość dynamiczna czynnika chłodniczego (ciecz i gaz)
μ_0l = CP.PropsSI('V','Q',0,'T',T_0,refrigerant)
print ("* Lepkość dynamiczna %s (ciecz): %.6f [Pas]" % (refrigerant,μ_0l))
μ_0v = CP.PropsSI('V','Q',1,'T',T_0,refrigerant)
print ("* Lepkość dynamiczna %s (gaz): %.6f [Pas]" % (refrigerant,μ_0v))

# współczynnik przewodzenia ciepła
λ_0l = CP.PropsSI('L','Q',0,'T',T_0,refrigerant)
print ("* Współczynnik przewodzenia ciepła (ciecz) %s: %.3f [W/mK]" % (refrigerant,λ_0l))
λ_0v = CP.PropsSI('L','Q',1,'T',T_0,refrigerant)
print ("* Współczynnik przewodzenia ciepła (gaz) %s: %.3f [W/mK]" % (refrigerant,λ_0v))

# napięcie powierzchniowe
σ_0 = CP.PropsSI('I','Q',0,'T',T_0,refrigerant)
print ("* Napięcie powierzchniowe ciepła %s: %.3f [N/m]" % (refrigerant,σ_0))

# Liczba Prandtla
Pr_0l = c_p0l * μ_0l / λ_0l
print ("* Liczba Prandtla (ciecz) %s: %.3f [-]" % (refrigerant,Pr_0l))
Pr_0v = c_p0v * μ_0v / λ_0v
print ("* Liczba Prandtla (gaz) %s: %.3f [-]" % (refrigerant,Pr_0l))


# ### Obliczanie współczynnika wnikania ciepła po stronie wrzącego czynnika
# Współczynnik wnikania ciepła podczas wrzenia w rurkach poziomych obliczony jest na podstawie metody Lavina i Younga.
# 
# **Źródło:** http://onlinelibrary.wiley.com/doi/10.1002/aic.690110629/abstract

# Obliczenie współczynnika wnikania ciepła po stronie czynnika jest możliwe przy znajomości (wcześniejszym obliczeniu/założeniu) następujących parametrów:
# * średnicy rurek $d_w$
# * strumienia masy czynnika chłodniczego $m_0$

# In[7]:


# Wymiary rurek, średnica zewnętrzna, wewnętrzna i grubość ścianki) oraz współczynnik przewodzenia ciepła materiału
# rurki
d_z = 0.01 # [m] - średnica zewnętrzna rurki
δ_r = 0.001 # [m] - grubość ścianki rurki
d_w = d_z - 2 * δ_r # [m] - średnica wewnętrzna rurki
λ_r = 370 # [W/mK] - miedź


# In[8]:


# ilość strumieni masy czynnika
n_strumieni = 2

# Frakcja ciekłego czynnika chłodniczego opuszczającego wymiennik
x = 0.5


# In[9]:


# Gęstość strumienia masy czynnika w rurkach
G_0 = (4 * m_0) / (n_strumieni * math.pi * d_w**2)
print ("* Gęstość strumienia masy czynnika: %.0f [kg/m2s]" % G_0)

# Liczba Reynoldsa obliczona dla frakcji przepływającej cieczy
Re_0 = G_0 * (1 - x) * d_w / μ_0l
print ("* Liczba Reynoldsa obliczona dla frakcji przepływającej cieczy: %.0f" % Re_0)


# Wartość liczby Nusselta okresla się wg korelacji Siedler-Tate'a. Należy pamiętać, ze korelacja ta znajduje potwierdzenie dla przepływów turbulentnych (Re > 10000).
# 
# $$Nu = 0.023 \cdot Re^{0.8} \cdot Pr^{1/3} \cdot \left( \frac{\mu_l}{\mu_{sc}} \right)^{0.14}$$
# gdzie:
# * $\mu_l$ - lepkość ciekłego czynnika w temperaturze wrzenia $T_0$
# * $\mu_{sc}$ - lepkość ciekłego czynnika w temperaturze ścianki, którą można przyjąć jako równą $T_0 + LMTD$
# 

# In[10]:


# lepkość płynu w temperaturze ścianki, tutaj przyjętej jako średnia logarytmiczna różnica temperatur LMTD
μ_sc =  CP.PropsSI('V','Q',0,'T',T_0 + LMTD,refrigerant)
print ("* Lepkość dynamiczna %s w temperaturze ścianki %.2f [K]: %.6f [Pas]" % (refrigerant,T_0+LMTD,μ_sc))

Nu_0 = 0.023 * Re_0**0.8 * Pr_0l**(1/3) * (μ_0l/μ_sc)**0.14
print ("* Liczba Nusselta: %.2f " % Nu_0)

# Współczynnik wnikania ciepła dla założonej frakcji cieczy
α_0l = Nu_0 * λ_0l / d_w
print ("* Współczynnik wnikania ciepła dla założonej frakcji cieczy: %.2f [W/m2K]" % α_0l)


# #### Obliczanie wartości współczynnika wnikania ciepła podczas przepływu mieszaniny cieczowo-parowej.
# Współczynnik wnikania ciepła dla mieszaniny cieczowo-parowej. Stała C w równaniu jest to współczynnik zależny od typu i orientacji rurki. Dla rur gładkich instalowanych poziomo C = 6.59, natomiast dla rur gładkich instalowanych pionowo C = 3.79.

# In[11]:


C = 6.59


# In[12]:


q = 1000 # wartość początkowa w iteracji
α_it = α_0l

q_old = 0
print (x)
while abs(q - q_old) > 1e-9:
    q_old = q
    print ("q = %.3f [W/m2]" % (q))    
    
    q = LMTD / ( (1/α_it) + (δ_r/λ_r) + (1/α_p))
    α_it = ( (C * ((1+x) / (1-x))**1.16) / ((G_0 * r_0) / q)**0.1 ) * α_0l

    print ("* Współczynnik przejmowania ciepła: %.3f [W/m2K]" % α_it)

α_0 = α_it


# In[ ]: