{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6de1cde8",
   "metadata": {},
   "source": [
    "# Heat Pump Model\n",
    "\n",
    "Here is a simple example of how to do a heat pump cycle calculation\n",
    "for a simple four-component system, with very simple models\n",
    "for each component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "46518868",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:10:59.064818Z",
     "iopub.status.busy": "2024-12-12T18:10:59.064658Z",
     "iopub.status.idle": "2024-12-12T18:10:59.353349Z",
     "shell.execute_reply": "2024-12-12T18:10:59.352834Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'name': 'R125',\n",
       " 'COP': np.float64(3.238624644077594),\n",
       " 'pevap / kPa': np.float64(606.2416455945954),\n",
       " 'pcond / kPa': np.float64(2001.2735289000334),\n",
       " 'rho1 / mol/m^3': np.float64(306.9219167132696),\n",
       " 'Qvol / kJ/m^3': np.float64(3287.7920883249403)}"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from dataclasses import dataclass \n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import teqp\n",
    "\n",
    "R = 8.31446261815324\n",
    "k_B = 1.380649e-23\n",
    "N_A = R/k_B\n",
    "\n",
    "@dataclass\n",
    "class ModelCombo:\n",
    "    \n",
    "    model : object\n",
    "    aig: object\n",
    "    name: str\n",
    "\n",
    "    def get_h(self, T, rhomolar, molefrac):\n",
    "        Atot10 = self.model.get_Ar10(T, rhomolar, molefrac) + self.aig.get_Aig10(T, rhomolar, molefrac)\n",
    "        return R*T*(1+self.model.get_Ar01(T, rhomolar, molefrac) + Atot10)\n",
    "\n",
    "    def get_s(self, T, rhomolar, molefrac):\n",
    "        Atot10 = self.model.get_Ar10(T, rhomolar, molefrac) + self.aig.get_Aig10(T, rhomolar, molefrac)\n",
    "        Atot00 = self.model.get_Ar00(T, rhomolar, molefrac) + self.aig.get_Aig00(T, rhomolar, molefrac)\n",
    "        return R*(Atot10 - Atot00)\n",
    "\n",
    "    def get_p(self, T, rhomolar, molefrac):\n",
    "        return rhomolar*self.model.get_R(molefrac)*T*(1 + self.model.get_Ar01(T, rhomolar, molefrac))\n",
    "\n",
    "def cycle(combo, anc, *, Tevap, Tcond, DELTAT_sh, DELTAT_sc, eta_comp):\n",
    "    \"\"\"\n",
    "    combo(ModelCombo): The joined model with ideal-gas and residual portions\n",
    "    anc : A set of ancillary functions implementing rhoL(T) and rhoV(T) methods\n",
    "    Tevap : Saturated vapor temperature in evaporator, in K\n",
    "    Tcond : Saturated vapor temperature in condenser, in K\n",
    "    DELTAT_sh : superheat, in K\n",
    "    DELTAT_sc : subcooling, in K\n",
    "    eta_comp : compressor efficiency\n",
    "    \"\"\"\n",
    "\n",
    "    model = combo.model\n",
    "    z = np.array([1.0])\n",
    "\n",
    "    # VLE densities,\n",
    "    # w/ guess values from the ancillary\n",
    "    rhomolar1satL, rhomolar1sat = model.pure_VLE_T(Tevap, anc.rhoL(Tevap), anc.rhoV(Tevap), 10)\n",
    "    rhomolar3sat, rhomolar3satV = model.pure_VLE_T(Tcond, anc.rhoL(Tcond), anc.rhoV(Tcond), 10)\n",
    "    p1 = combo.get_p(Tevap, rhomolar1sat, z)\n",
    "    p2 = combo.get_p(Tcond, rhomolar3sat, z)\n",
    "\n",
    "    # Evaporator outlet & compressor inlet @ state point 1\n",
    "    T1 = Tevap + DELTAT_sh\n",
    "    rhomolar1 = scipy.optimize.newton(lambda rho_: combo.get_p(T1, rho_, z)-p1, rhomolar1sat)\n",
    "    h1 = combo.get_h(T1, rhomolar1, z)\n",
    "    s1 = combo.get_s(T1, rhomolar1, z)\n",
    "    \n",
    "    # Solve for isentropic compressor outlet\n",
    "    res = lambda x: [combo.get_p(x[0], x[1], z)-p2, combo.get_s(x[0], x[1], z)-s1]\n",
    "    T2s, rho2s = scipy.optimize.fsolve(res, [Tcond, rhomolar1sat])\n",
    "    h2s = combo.get_h(T2s, rho2s, z)\n",
    "    h2 = h1 + (h2s-h1)/eta_comp # @ state point 2\n",
    "\n",
    "    # Condenser outlet and expansion valve inlet @ state point 3\n",
    "    T3 = Tcond - DELTAT_sc\n",
    "    rhomolar3 = scipy.optimize.newton(lambda rho_: combo.get_p(T3, rho_, z)-p2, rhomolar3sat)\n",
    "    h3 = combo.get_h(T3, rhomolar3, z)\n",
    "\n",
    "    COP = (h1-h3)/(h2-h1)\n",
    "\n",
    "    return {\n",
    "        'name': combo.name,\n",
    "        'COP': COP,\n",
    "        'pevap / kPa': p1/1e3, \n",
    "        'pcond / kPa': p2/1e3,\n",
    "        'rho1 / mol/m^3': rhomolar1,\n",
    "        'Qvol / kJ/m^3': (h1-h3)*rhomolar1/1e3,\n",
    "    }\n",
    "\n",
    "# Build the model (ideal-gas and residual)\n",
    "FLD = 'R125'\n",
    "path = teqp.get_datapath()+f'/dev/fluids/{FLD}.json'\n",
    "assert(os.path.exists(path))\n",
    "jig = teqp.convert_CoolProp_idealgas(path, 0)\n",
    "combo = ModelCombo(\n",
    "    model=teqp.build_multifluid_model([path], ''),\n",
    "    aig=teqp.IdealHelmholtz([jig]),\n",
    "    name=FLD\n",
    ")\n",
    "\n",
    "# Generic ancillary functions from on-the-fly ancillary construction\n",
    "anc = teqp.build_ancillaries(combo.model, 360, 6000, 250)\n",
    "\n",
    "cycle(combo, anc=anc, Tevap=270, Tcond=313, DELTAT_sh=5, DELTAT_sc=5, eta_comp=0.7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4be7c14a",
   "metadata": {},
   "source": [
    "Exercise for the reader: plot the points on a P-H diagram, showing the saturated liquid and vapor states"
   ]
  }
 ],
 "metadata": {
  "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}