{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f225380f",
   "metadata": {},
   "source": [
    "# Example: Departure functions from Peng-Robinson equation of state\n",
    "\n",
    "These are codes to calculate the enthalpy departure function from the complete, generalized Peng-Robinson equation of state as given in SIS 5th edition equations 6.4-2 and 6.7-1 to 6.7-4. \n",
    "\n",
    "We use it to calculate analyze an isenthalpic expansion of methane through a throttling process.\n",
    "\n",
    "Eric Furst  \n",
    "November 2024\n",
    "\n",
    "Version updates:  \n",
    "Revisions in layout, calculation of PR EOS (11/2025)\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28594ad2",
   "metadata": {},
   "source": [
    "As ususal, we import the *numpy* and *matplotlib* libraries. We'll also use the *SciPy* constants library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8ba8e9f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import constants\n",
    "import scipy.optimize as opt\n",
    "\n",
    "R = constants.R # Set the gas constant to R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8672870c-eb49-47d7-b862-437a5b5515d6",
   "metadata": {},
   "source": [
    "## Problem statement\n",
    "\n",
    "Methane gas undergoes a rapid, adiabatic expansion through a continous throttling process from upstream conditions ${\\text{13}^\\circ}\\text{C}$ and 184 bar to a lower pressure at ${\\text{-43}^\\circ}\\text{C}$ .\n",
    "\n",
    "Calculate the downstream gas pressure.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2e66192-d13c-4c0d-9677-6a0e88410ff3",
   "metadata": {},
   "source": [
    "To solve this problem, we will calculate the enthalpy of methane such that it satisfies an isenthalpic condition. In terms of the departure function values, this is \n",
    "$$\\Delta \\underline{H} = \\Delta \\underline{H}^\\mathrm{IG} + \\Delta [ \\underline{H} - \\underline{H}^\\mathrm{IG}] = 0$$\n",
    "or\n",
    "$$\\Delta \\underline{H}^\\mathrm{IG} + [ \\underline{H} - \\underline{H}^\\mathrm{IG}]_2 - [ \\underline{H} - \\underline{H}^\\mathrm{IG}]_1 = 0$$\n",
    "The specific departure functions we will use are derived from the generalized Peng-Robinson equation of state. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cacacf5",
   "metadata": {},
   "source": [
    "# Generalized Peng-Robinson Equation of State\n",
    "\n",
    "The Generalized Peng-Robinson equation of state is\n",
    "\n",
    "$$ P = \\frac{RT}{\\underline{V}-b} - \\frac{a(T)}{\\underline{V}(\\underline{V}+b) + b(\\underline{V}-b)} \\tag{Eq. 6.4-2}$$\n",
    "\n",
    "with\n",
    "\n",
    "$$ b = 0.07780 \\frac{RT_c}{P_c} \\tag{Eq. 6.7-2}$$\n",
    "$$ a(T) = a(T_c)\\alpha(T) = 0.45724 \\frac{R^2T_c^2}{P_c}\\alpha(T) \\tag{Eq. 6.7-1}$$\n",
    "$$ \\sqrt{\\alpha} = 1 + \\kappa \\left ( 1- \\sqrt{\\frac{T}{T_c}} \\right ) \\tag{Eq. 6.7-3}$$\n",
    "$$ \\kappa = 0.37464 + 1.54226\\omega − 0.26992\\omega^2 \\tag{Eq. 6.7-4}$$\n",
    "\n",
    "The acentric factor $\\omega$ and the crticial temperatures and pressures are given in SIS table 6.6-1.\n",
    "\n",
    "Calculating the pressure $P$ given $\\underline{V}$ and $T$ is straightforward, but to calculate the molar volume given $P$ and $T$, we need to solve the cubic equation of state of the form\n",
    "\n",
    "$$ Z^3 + \\alpha Z^2 + \\beta Z + \\gamma = 0 \\tag{Eq. 6.4-4}$$\n",
    "\n",
    "where $Z$ is the compressibility factor\n",
    "\n",
    "$$ Z = \\frac{P \\underbar{V{}}}{RT} $$\n",
    "\n",
    "For the Peng-Robinson EOS,\n",
    "\n",
    "$$ \\alpha = -1 + B $$\n",
    "$$ \\beta = A - 3B^2 -2B $$\n",
    "$$ \\gamma = -AB + B^2 + B^3 $$\n",
    "\n",
    "and \n",
    "\n",
    "$$ A = \\frac{aP}{(RT)^2} $$\n",
    "$$ B = \\frac{bP}{RT} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f3d2f917",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Generalized Peng-Robinson EOS \n",
    "PR_pressure returns the pressure given V, T, Pc, Tc, omega\n",
    "PR_volume returns all molar volumes (real roots of EOS) given P, T, Pc, Tc, omega\n",
    "\"\"\"\n",
    "\n",
    "from scipy import constants\n",
    "from numpy.polynomial import Polynomial\n",
    "\n",
    "R = constants.R # Set the gas constant to R\n",
    "\n",
    "def calc_b(Pc,Tc):\n",
    "    return 0.07780*R*Tc/Pc\n",
    "\n",
    "def calc_a(T,Pc,Tc,omega):\n",
    "    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2\n",
    "    sqrtalpha = 1 + kappa*(1-np.sqrt(T/Tc))\n",
    "    return 0.45724*R**2*Tc**2/Pc*sqrtalpha**2\n",
    "\n",
    "# Calculate the presure given V, T for PR EOS\n",
    "def PR_pressure(V,T,Pc,Tc,omega):\n",
    "    a = calc_a(T,Pc,Tc,omega)\n",
    "    b = calc_b(Pc,Tc)\n",
    "\n",
    "    P = R*T/(V-b) - a/(V*(V+b)+b*(V-b))\n",
    "    return P\n",
    "\n",
    "# Solve for the compressibility factor given P, T for PR EOS in m^3/mol\n",
    "# Note that we can return multiple real roots (up to three)\n",
    "# The largest and smallest will be the vapor and liquid, respectively\n",
    "def PR_compressibility(P, T, Pc, Tc, omega):\n",
    "    # Calculate a, b, A, and B\n",
    "    a = calc_a(T, Pc, Tc, omega)\n",
    "    b = calc_b(Pc, Tc)\n",
    "    A = a*P/R**2/T**2\n",
    "    B = b*P/R/T\n",
    "\n",
    "    # Definitions of alpha, beta, gamma in SIS Table 6.4-3 for PR EOS\n",
    "    alpha = -1 + B\n",
    "    beta = A - 3*B**2 -2*B\n",
    "    gamma = -A*B + B**2 + B**3\n",
    "\n",
    "    # polynomial with coefficients in increasing order: c0 + c1 x + c2 x**2 + ...\n",
    "    p = Polynomial([ gamma, beta, alpha, 1 ])  \n",
    "\n",
    "    roots = p.roots()        # returns all (possibly complex) roots of Z\n",
    "    real_roots = roots.real[abs(roots.imag) < 1e-12]  # filter real ones\n",
    "\n",
    "    return real_roots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74031b60-3999-49f8-91a2-eebfeeb09d83",
   "metadata": {},
   "source": [
    "## Data for methane\n",
    "\n",
    "For methane, we need the critical pressure and temperature and Pitzer's acentric factor:\n",
    "\n",
    "$P_c = 4.6$ MPa   \n",
    "$T_c = 190.6$ K  \n",
    "$\\underline{V}_c = 0.099\\,\\mathrm{m^3/kmol} = 9.9\\times10^{-5}\\,\\mathrm{m^3/mol}$  \n",
    "$\\omega = 0.008$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7579e5f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Critical values for Methane and the accentric factor\n",
    "# pressure in Pa, temperature in K\n",
    "Pc = 4.6e6\n",
    "Tc = 190.6\n",
    "omega = 0.008\n",
    "\n",
    "# Known values in problem\n",
    "T1 = 286        # Inlet temperature in K\n",
    "P1 = 18.4e6     # Inlet pressure in Pa\n",
    "T2 = 230        # Outlet temperature in K"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad238854-9629-4ce9-9f2a-0c92edcecc92",
   "metadata": {},
   "source": [
    "# Enthalpy calculation\n",
    "\n",
    "## Start by calculating $Z$ for the inlet condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0bdc3b7d-bbe1-4358-91a2-8aeaa0d33040",
   "metadata": {},
   "outputs": [],
   "source": [
    "TR1 = T1/Tc      # Reduced temperature\n",
    "PR1 = P1/Pc      # Reduced pressure\n",
    "\n",
    "# Use the PR_vol function to calculate the molar volume\n",
    "Z1 = np.max(PR_compressibility(P1,T1,Pc,Tc,omega))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ca49e458-2711-465b-add8-198d80ff5131",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z1 = 0.77\n"
     ]
    }
   ],
   "source": [
    "print(\"Z1 = {:.2f}\".format(Z1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4133d8cd",
   "metadata": {},
   "source": [
    "## Calculate the departure function for the inlet condition\n",
    "The enthalpy departure function in $(P,T)$ control variables for the Peng-Robinson equuation of state is\n",
    "\n",
    "$$[\\underline{H} - \\underline{H}^\\mathrm{IG}] = RT(Z-1) + \\frac{T\\frac{da}{dT}-a}{2\\sqrt2 b}\\ln \\left [ \\frac{Z+(1+\\sqrt2)B}{Z+(1-\\sqrt2)B} \\right ] \\tag{Eq. 6.4-29}$$\n",
    "\n",
    "where we need to evaluate the compressibility factor,\n",
    "$$Z = \\frac{P\\underline{V}}{RT}$$\n",
    "the scaled value of $b$,\n",
    "$$B = \\frac{bP}{RT}$$\n",
    "and the deriviative of $a(T)$ with respect to temperature given in SIS (p. 264),\n",
    "$$\\frac{da}{dT} = - 0.45724 \\frac{R^2T_c^2}{P_c} \\kappa \\sqrt{\\frac{\\alpha}{T T_c}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b9c7dd01-8aca-45eb-9d06-7eb0d1d1c084",
   "metadata": {},
   "outputs": [],
   "source": [
    "# With Z already calculated, we define a function to calculate the departure function\n",
    "\n",
    "\"\"\"\n",
    "Enthalpy departure function from the Peng-Robinson equation of state (eq. 6.4-29)\n",
    "For convenience, we recalculate b, kappa, sqrtalpha, and a(T) here, but we could\n",
    "use the earlier functions, too.\n",
    "\n",
    "Uses globals Tc, Pc, R\n",
    "\"\"\"\n",
    "def calc_depH(T,P,Z):\n",
    "    b = 0.07780*R*Tc/Pc\n",
    "    B = b*P/R/T\n",
    "    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2\n",
    "    sqrtalpha = 1 + kappa*(1-(T/Tc)**0.5)\n",
    "    a = 0.45724*R**2*Tc**2/Pc*sqrtalpha**2\n",
    "    dadT = -0.45724*R*R*Tc*Tc/Pc*kappa*sqrtalpha/(T*Tc)**0.5\n",
    "\n",
    "    depH = R*T*(Z-1)+(T*dadT-a)/(2*2**0.5*b)*np.log((Z+(1+2**0.5)*B)/(Z+(1-2**0.5)*B))\n",
    "    return depH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c3dc8ef7-f035-4079-8e90-bb4cc500fbf9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** depH1 = -3134 J/mol\n"
     ]
    }
   ],
   "source": [
    "depH1 = calc_depH(T1,P1,Z1)\n",
    "\n",
    "print(f\"*** depH1 = {depH1:.0f} J/mol\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6132bf00-dc64-4bbf-99cc-4e749d5f676b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.93\n"
     ]
    }
   ],
   "source": [
    "# I'm calculating this number to compare with corresponding states later...\n",
    "print(f\"{-depH1/Tc/4.184:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b0c2408-14d8-4e9b-b593-414804278b6d",
   "metadata": {},
   "source": [
    "## Calculate the change in enthalpy in the ideal state\n",
    "$$\\Delta\\underline{H}^\\mathrm{IG} = \\int_{T_1}^{T_2}C_p^*(T)dT$$\n",
    "where we will use the ideal state heat capacities summarized in Appendix A.II,\n",
    "$$C^*_p(T) = A + BT + CT^2 + DT^3$$\n",
    "resulting in\n",
    "$$\\Delta\\underline{H}^\\mathrm{IG} = A(T_2-T_1) + \\frac{1}{2}B(T_2^2 - T_1^2)+ \\frac{1}{3}C(T_2^3 - T_1^3) + \\frac{1}{4}D(T_2^4 - T_1^4)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b587ff94-888c-455a-b216-f68a3a7ff462",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heat capacity data for N2, O2, CO2, CH4, C2H6\n",
    "params = np.array([\n",
    "    [28.83, -0.157, 0.808, -2.871, 1800],\n",
    "    [25.46, 1.519, -0.715, 1.311, 1800],\n",
    "    [22.243, 5.977, -3.499, 7.464, 1800],\n",
    "    [19.875, 5.021, 1.268, -11.004, 1500],\n",
    "    [6.895, 17.255, -6.402, 7.280, 1500]\n",
    "    ])\n",
    "\n",
    "i = params[3,:]  # equation parameters for methane\n",
    "\n",
    "deltaHIG = i[0]*(T2-T1) + i[1]*10**-2*(T2**2-T1**2)/2 \\\n",
    "            + i[2]*10**-5*(T2**3-T1**3)/3 \\\n",
    "            + i[3]*10**-9*(T2**4-T1**4)/4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bb8d3194-40ce-463c-9914-5eed37d4f7de",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Delta HIG = -1875 J/mol\n"
     ]
    }
   ],
   "source": [
    "print(\"Delta HIG = {:.0f} J/mol\".format(deltaHIG))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d72851d-90a5-43f0-9839-02f5f495744b",
   "metadata": {},
   "source": [
    "Now solve\n",
    "$$[ \\underline{H} - \\underline{H}^\\mathrm{IG}]_2 = [ \\underline{H} - \\underline{H}^\\mathrm{IG}]_1 - \\Delta \\underline{H}^\\mathrm{IG}$$\n",
    "\n",
    "We are looking for a value of the departure function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9907b191-2e60-4a7b-8f04-e2655ed4474c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "depH2 must equal -1259 J/mol\n"
     ]
    }
   ],
   "source": [
    "print(\"depH2 must equal {:.0f} J/mol\".format(depH1-deltaHIG))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5cf19a3-a21e-4fc0-a034-30d48834da91",
   "metadata": {},
   "source": [
    "Do this by guessing an outlet pressure, calculating $[ \\underline{H} - \\underline{H}^\\mathrm{IG}]_2$ and iterating (guessing new values of $P_2$) until the condition above is met."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04bda783-43d6-4fb6-8f69-e87dc1a5a2be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "depH2 = -1259 J/mol, deltaH = 0 J/mol <-- should be zero\n"
     ]
    }
   ],
   "source": [
    "# Guess:\n",
    "P2 = 4.145e6  # guess a pressure\n",
    "\n",
    "# Use the PR_compressibility function to calculate the molar volume and compressibility factor\n",
    "Z2 = np.max(PR_compressibility(P2,T2,Pc,Tc,omega))\n",
    "\n",
    "depH2 = calc_depH(T2,P2,Z2)\n",
    "#print(P2/Pc)\n",
    "print(\"depH2 = {:.0f} J/mol, deltaH = {:.0f} J/mol <-- should be zero\".format(depH2,deltaHIG+depH2-depH1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acfea36c-cb4a-4a12-8577-059d4790340e",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Change in molar entropy\n",
    "\n",
    "Now we can calcuate the change in molar entropy of the methane,\n",
    "\n",
    "$$\\Delta \\underline{S} = \\Delta \\underline{S}^\\mathrm{IG} + \\Delta [ \\underline{S} - \\underline{S}^\\mathrm{IG}]$$\n",
    "\n",
    "Since the throttling process is adiabatic, this change represents the rate of entropy generation, $\\underline{\\dot{S}}_\\mathrm{gen}$.\n",
    "\n",
    "The entropy departure function for the Peng-Robinson equation is\n",
    "$$[ \\underline{S} - \\underline{S}^\\mathrm{IG}] = R\\ln(Z-B) + \\frac{da/dT}{2\\sqrt{2}b}\\ln \\left [ \\frac{Z+(1+\\sqrt2)B}{Z+(1-\\sqrt2)B} \\right ] \\tag{Eq. 6.4-30}$$\n",
    "and the entropy change of the real fluid in an ideal state with $P$ and $T$ control variables is\n",
    "$$\\Delta \\underline{S}^\\mathrm{IG} = \\int_{T_1}^{T_2}\\frac{C_p^*(T)}{T}dT - R\\ln\\frac{P_2}{P_1}$$\n",
    "with $C_p^*(T)$ given above, \n",
    "$$\\Delta \\underline{S}^\\mathrm{IG} = A\\ln\\frac{T_2}{T_1} + B(T_2-T_1) + \\frac{1}{2}C(T_2^2-T_1^2) + \\frac{1}{3}D(T_2^3-T_1^3) - R\\ln\\frac{P_2}{P_1}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cea15f29-14c8-4b4e-91bc-b070b56856d1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(12.433452523746052)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "deltaSIG = i[0]*np.log(T2/T1) + i[1]*10**-2*(T2-T1) + i[2]*10**-5*(T2**2-T1**2)/2 \n",
    "+ i[3]*10**-9*(T2**3-T1**3)/3 - R*np.log(P2/P1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0050af8f-231e-401a-8583-4e536511082a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Entropy departure function from the Peng-Robinson equation of state (eq. 6.4-30)\n",
    "For convenience, we recalculate b, kappa, sqrtalpha, and a(T) here, but we could\n",
    "use the earlier functions, too.\n",
    "\n",
    "Uses globals Tc, Pc, R\n",
    "\"\"\"\n",
    "def calc_depS(T,P,Z):\n",
    "    b = 0.07780*R*Tc/Pc\n",
    "    B = b*P/R/T\n",
    "    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2\n",
    "    sqrtalpha = 1 + kappa*(1-np.sqrt(T/Tc))\n",
    "    a = 0.45724*R**2*Tc**2/Pc*sqrtalpha**2\n",
    "    dadT = -0.45724*R*R*Tc*Tc/Pc*kappa*sqrtalpha/np.sqrt(T*Tc)\n",
    "\n",
    "    depS = R*(Z-B)+dadT/(2*2**0.5*b)*np.log((Z+(1+2**0.5)*B)/(Z+(1-2**0.5)*B))\n",
    "    return depS\n",
    "\n",
    "depS1 = calc_depS(T1, P1, Z1)\n",
    "depS2 = calc_depS(T2, P2, Z2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f21e77c3-ac10-4596-bcd1-1e1124cffdd6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T1: 286.0 K, P1: 1.84e+07 Pa, Z1: 0.77\n",
      "T2: 230.0 K, P2: 4.14e+06 Pa, Z2: 0.79\n",
      "\n",
      "delSIG= -7.33, depS2= 4.9, depS1= 1.6\n",
      "\n",
      "Change in molar entropy: -4.0 J/mol.K\n"
     ]
    }
   ],
   "source": [
    "print(\"T1: {:.1f} K, P1: {:0.2e} Pa, Z1: {:.2f}\".format(T1,P1,Z1))\n",
    "print(\"T2: {:.1f} K, P2: {:0.2e} Pa, Z2: {:.2f}\\n\".format(T2,P2,Z2))\n",
    "print(\"delSIG= {:.2f}, depS2= {:.1f}, depS1= {:.1f}\\n\".format(deltaSIG,depS2,depS1))\n",
    "print(\"Change in molar entropy: {:.1f} J/mol.K\".format(deltaSIG + depS2 - depS1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2f846c5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}