{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "import CoolProp\n",
    "import CoolProp.CoolProp as CP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorytm obliczeń współczynnika wnikania ciepła po stronie wrzącego czynnika\n",
    "Podczas wrzenia czynnika w rurkach poziomych."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Założenia projektowe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_k = 15 * 1e3 # [kW] - wydajność grzewcza projektowanej pompy ciepła\n",
    "\n",
    "T_0 = 3 + 273.15 # [K] - temperatura wrzenia czynnika chłodniczego\n",
    "T_k = 40 + 273.15 # [K] - temperatura skraplania czynnika chłodniczego\n",
    "dT_sup = 10 # [°C] - przegrzanie\n",
    "\n",
    "n_iz = 0.5 # założony współczynnik sprawności izentropowej sprężarki"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# czynnik chłodniczy\n",
    "refrigerant = \"R134a\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Obliczanie punktów charakterystycznych obiegu pompy ciepła"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "* Strumień masy czynnika chłodniczego: 0.074 [kg/s]\n",
      "* Wydajność chłodnicza parownika: 11.312 [kW]\n"
     ]
    }
   ],
   "source": [
    "p_0 = CP.PropsSI('P','Q',1,'T',T_0,refrigerant)\n",
    "p_k = CP.PropsSI('P','Q',1,'T',T_k,refrigerant)\n",
    "\n",
    "if dT_sup > 0: h1 = CP.PropsSI('H','P',p_0,'T',T_0+dT_sup,refrigerant)\n",
    "else: h1 = CP.PropsSI('H','Q',1,'T',T_0,refrigerant)\n",
    "s1 = CP.PropsSI('S','H',h1,'P',p_0,refrigerant)\n",
    "h2 = CP.PropsSI('H','S',s1,'P',p_k,refrigerant)\n",
    "h3 = CP.PropsSI('H','Q',0,'T',T_k,refrigerant)\n",
    "h4 = h3\n",
    "\n",
    "T_2 = CP.PropsSI('T','S',s1,'P',p_k,refrigerant)\n",
    "\n",
    "# Korekta położenia punktu 2 po uwzględnieniu sprawności izentropowej sprężarki\n",
    "h2b = h1 + (h2 - h1) / n_iz\n",
    "\n",
    "q_0 = h1 - h4\n",
    "q_k = h2b - h3\n",
    "\n",
    "# strumień masy czynnika chłodniczego\n",
    "m_0 = Q_k / q_k\n",
    "print (\"* Strumień masy czynnika chłodniczego: %.3f [kg/s]\" % m_0)\n",
    "\n",
    "Q_0 = m_0 * q_0\n",
    "print (\"* Wydajność chłodnicza parownika: %.3f [kW]\" % (Q_0/1e3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dane podstawowe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_p1 = 12 + 273.15 # [K] - temperatura powietrza na wlocie do parownika\n",
    "T_p2 = 7 + 273.15 # [K] - temperatura powietrza na wylocie z parownika\n",
    "\n",
    "# Średnia logarytmiczna różnica temperatur wymiennika\n",
    "LMTD = (T_p1 - T_0) - (T_p2 - T_0) / math.log((T_p1 - T_0)/(T_p2 - T_0))\n",
    "\n",
    "# Współczynnik wnikania ciepła po stronie powietrza (powinien być wcześniej obliczony)\n",
    "α_p = 70 # [W/m2K]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Czynnik chłodniczy i jego parametry w założonych temperaturach pracy wymiennika"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "* Ciśnienie nasycenia R134a: 325984.93 [J/kgK]\n",
      "* Ciepło właściwe R134a (ciecz): 1349.39 [J/kgK]\n",
      "* Ciepło właściwe R134a (gaz): 911.08 [J/kgK]\n",
      "* Ciepło przemiany fazowej R134a: 196301.10 [J/kg]\n",
      "* Gęstość R134a (ciecz): 1284.80 [kg/m3]\n",
      "* Gęstość R134a (gaz): 16.01 [kg/m3]\n",
      "* Lepkość dynamiczna R134a (ciecz): 0.000257 [Pas]\n",
      "* Lepkość dynamiczna R134a (gaz): 0.000011 [Pas]\n",
      "* Współczynnik przewodzenia ciepła (ciecz) R134a: 0.091 [W/mK]\n",
      "* Współczynnik przewodzenia ciepła (gaz) R134a: 0.012 [W/mK]\n",
      "* Napięcie powierzchniowe ciepła R134a: 0.011 [N/m]\n",
      "* Liczba Prandtla (ciecz) R134a: 3.817 [-]\n",
      "* Liczba Prandtla (gaz) R134a: 3.817 [-]\n"
     ]
    }
   ],
   "source": [
    "# ciśnienie nasycenia czynnika chłodniczego\n",
    "p_n0 = CP.PropsSI('P','Q',0,'T',T_0,refrigerant)\n",
    "print (\"* Ciśnienie nasycenia %s: %.2f [J/kgK]\" % (refrigerant,p_n0))\n",
    "\n",
    "# ciepło właściwe czynnika chłodniczego\n",
    "c_p0l = CP.PropsSI('C','Q',0,'T',T_0,refrigerant)\n",
    "print (\"* Ciepło właściwe %s (ciecz): %.2f [J/kgK]\" % (refrigerant,c_p0l))\n",
    "c_p0v = CP.PropsSI('C','Q',1,'T',T_0,refrigerant)\n",
    "print (\"* Ciepło właściwe %s (gaz): %.2f [J/kgK]\" % (refrigerant,c_p0v))\n",
    "\n",
    "\n",
    "# ciepło przemiany fazowej\n",
    "r_0 = CP.PropsSI('H','Q',1.0,'T',T_0,refrigerant) - CP.PropsSI('H','Q',0,'T',T_0,refrigerant)\n",
    "print (\"* Ciepło przemiany fazowej %s: %.2f [J/kg]\" % (refrigerant,r_0))\n",
    "\n",
    "# gęstość czynnika chłodniczego (ciecz i gaz)\n",
    "ρ_0l = CP.PropsSI('D','Q',0,'T',T_0,refrigerant)\n",
    "print (\"* Gęstość %s (ciecz): %.2f [kg/m3]\" % (refrigerant,ρ_0l))\n",
    "ρ_0v = CP.PropsSI('D','Q',1,'T',T_0,refrigerant)\n",
    "print (\"* Gęstość %s (gaz): %.2f [kg/m3]\" % (refrigerant,ρ_0v))\n",
    "\n",
    "# lepkość dynamiczna czynnika chłodniczego (ciecz i gaz)\n",
    "μ_0l = CP.PropsSI('V','Q',0,'T',T_0,refrigerant)\n",
    "print (\"* Lepkość dynamiczna %s (ciecz): %.6f [Pas]\" % (refrigerant,μ_0l))\n",
    "μ_0v = CP.PropsSI('V','Q',1,'T',T_0,refrigerant)\n",
    "print (\"* Lepkość dynamiczna %s (gaz): %.6f [Pas]\" % (refrigerant,μ_0v))\n",
    "\n",
    "# współczynnik przewodzenia ciepła\n",
    "λ_0l = CP.PropsSI('L','Q',0,'T',T_0,refrigerant)\n",
    "print (\"* Współczynnik przewodzenia ciepła (ciecz) %s: %.3f [W/mK]\" % (refrigerant,λ_0l))\n",
    "λ_0v = CP.PropsSI('L','Q',1,'T',T_0,refrigerant)\n",
    "print (\"* Współczynnik przewodzenia ciepła (gaz) %s: %.3f [W/mK]\" % (refrigerant,λ_0v))\n",
    "\n",
    "# napięcie powierzchniowe\n",
    "σ_0 = CP.PropsSI('I','Q',0,'T',T_0,refrigerant)\n",
    "print (\"* Napięcie powierzchniowe ciepła %s: %.3f [N/m]\" % (refrigerant,σ_0))\n",
    "\n",
    "# Liczba Prandtla\n",
    "Pr_0l = c_p0l * μ_0l / λ_0l\n",
    "print (\"* Liczba Prandtla (ciecz) %s: %.3f [-]\" % (refrigerant,Pr_0l))\n",
    "Pr_0v = c_p0v * μ_0v / λ_0v\n",
    "print (\"* Liczba Prandtla (gaz) %s: %.3f [-]\" % (refrigerant,Pr_0l))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Obliczanie współczynnika wnikania ciepła po stronie wrzącego czynnika\n",
    "Współczynnik wnikania ciepła podczas wrzenia w rurkach poziomych obliczony jest na podstawie metody Lavina i Younga.\n",
    "\n",
    "**Źródło:** http://onlinelibrary.wiley.com/doi/10.1002/aic.690110629/abstract"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:\n",
    "* średnicy rurek $d_w$\n",
    "* strumienia masy czynnika chłodniczego $m_0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wymiary rurek, średnica zewnętrzna, wewnętrzna i grubość ścianki) oraz współczynnik przewodzenia ciepła materiału\n",
    "# rurki\n",
    "d_z = 0.01 # [m] - średnica zewnętrzna rurki\n",
    "δ_r = 0.001 # [m] - grubość ścianki rurki\n",
    "d_w = d_z - 2 * δ_r # [m] - średnica wewnętrzna rurki\n",
    "λ_r = 370 # [W/mK] - miedź"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ilość strumieni masy czynnika\n",
    "n_strumieni = 2\n",
    "\n",
    "# Frakcja ciekłego czynnika chłodniczego opuszczającego wymiennik\n",
    "x = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "* Gęstość strumienia masy czynnika: 736 [kg/m2s]\n",
      "* Liczba Reynoldsa obliczona dla frakcji przepływającej cieczy: 11469\n"
     ]
    }
   ],
   "source": [
    "# Gęstość strumienia masy czynnika w rurkach\n",
    "G_0 = (4 * m_0) / (n_strumieni * math.pi * d_w**2)\n",
    "print (\"* Gęstość strumienia masy czynnika: %.0f [kg/m2s]\" % G_0)\n",
    "\n",
    "# Liczba Reynoldsa obliczona dla frakcji przepływającej cieczy\n",
    "Re_0 = G_0 * (1 - x) * d_w / μ_0l\n",
    "print (\"* Liczba Reynoldsa obliczona dla frakcji przepływającej cieczy: %.0f\" % Re_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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).\n",
    "\n",
    "$$Nu = 0.023 \\cdot Re^{0.8} \\cdot Pr^{1/3} \\cdot \\left( \\frac{\\mu_l}{\\mu_{sc}} \\right)^{0.14}$$\n",
    "gdzie:\n",
    "* $\\mu_l$ - lepkość ciekłego czynnika w temperaturze wrzenia $T_0$\n",
    "* $\\mu_{sc}$ - lepkość ciekłego czynnika w temperaturze ścianki, którą można przyjąć jako równą $T_0 + LMTD$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "* Lepkość dynamiczna R134a w temperaturze ścianki 280.22 [K]: 0.000244 [Pas]\n",
      "* Liczba Nusselta: 64.03 \n",
      "* Współczynnik wnikania ciepła dla założonej frakcji cieczy: 725.86 [W/m2K]\n"
     ]
    }
   ],
   "source": [
    "# lepkość płynu w temperaturze ścianki, tutaj przyjętej jako średnia logarytmiczna różnica temperatur LMTD\n",
    "μ_sc =  CP.PropsSI('V','Q',0,'T',T_0 + LMTD,refrigerant)\n",
    "print (\"* Lepkość dynamiczna %s w temperaturze ścianki %.2f [K]: %.6f [Pas]\" % (refrigerant,T_0+LMTD,μ_sc))\n",
    "\n",
    "Nu_0 = 0.023 * Re_0**0.8 * Pr_0l**(1/3) * (μ_0l/μ_sc)**0.14\n",
    "print (\"* Liczba Nusselta: %.2f \" % Nu_0)\n",
    "\n",
    "# Współczynnik wnikania ciepła dla założonej frakcji cieczy\n",
    "α_0l = Nu_0 * λ_0l / d_w\n",
    "print (\"* Współczynnik wnikania ciepła dla założonej frakcji cieczy: %.2f [W/m2K]\" % α_0l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Obliczanie wartości współczynnika wnikania ciepła podczas przepływu mieszaniny cieczowo-parowej.\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = 6.59"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5\n",
      "q = 1000.000 [W/m2]\n",
      "* Współczynnik przejmowania ciepła: 4557.009 [W/m2K]\n",
      "q = 259.630 [W/m2]\n",
      "* Współczynnik przejmowania ciepła: 4592.145 [W/m2K]\n",
      "q = 280.358 [W/m2]\n",
      "* Współczynnik przejmowania ciepła: 4592.199 [W/m2K]\n",
      "q = 280.390 [W/m2]\n",
      "* Współczynnik przejmowania ciepła: 4592.199 [W/m2K]\n",
      "q = 280.390 [W/m2]\n",
      "* Współczynnik przejmowania ciepła: 4592.199 [W/m2K]\n",
      "q = 280.390 [W/m2]\n",
      "* Współczynnik przejmowania ciepła: 4592.199 [W/m2K]\n"
     ]
    }
   ],
   "source": [
    "q = 1000 # wartość początkowa w iteracji\n",
    "α_it = α_0l\n",
    "\n",
    "q_old = 0\n",
    "print (x)\n",
    "while abs(q - q_old) > 1e-9:\n",
    "    q_old = q\n",
    "    print (\"q = %.3f [W/m2]\" % (q))    \n",
    "    \n",
    "    q = LMTD / ( (1/α_it) + (δ_r/λ_r) + (1/α_p))\n",
    "    α_it = ( (C * ((1+x) / (1-x))**1.16) / ((G_0 * r_0) / q)**0.1 ) * α_0l\n",
    "\n",
    "    print (\"* Współczynnik przejmowania ciepła: %.3f [W/m2K]\" % α_it)\n",
    "\n",
    "α_0 = α_it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}