{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f225380f",
   "metadata": {},
   "source": [
    "# Homework: Calculate the heat capacity of ethane at different pressures\n",
    "\n",
    "The argument we made for using departure functions is that the heat capacity for a real substance depends on pressure (or molar volume). Here we will calculate the constant pressure heat capacity at different pressures for ethane.\n",
    "\n",
    "SIS is Stanley I. Sandler, *Chemical, Biochemcial and Engineering Thermodynamics*, 5th ed.\n",
    "\n",
    "Eric Furst  \n",
    "November 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": 61,
   "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",
    "Using a Jupyter Notebook, plot the constant pressure heat capacity of ethane between 30°C and 200°C for the pressures 0.1, 2, and 4 MPa. Compare these to the ideal gas values $ C_p (T) $ given by the correlation in Appendix A.II."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2e66192-d13c-4c0d-9677-6a0e88410ff3",
   "metadata": {},
   "source": [
    "To solve this problem, I will calculate the enthalpy of ethane as a function of temperature at constant pressure. From these values, I will take the numerical derivative to find\n",
    "\n",
    "$$ C_p(T) = \\left ( \\frac{\\partial \\underbar{H}}{\\partial T} \\right )_P \\approx \\left ( \\frac{\\Delta \\underbar{H}}{\\Delta T} \\right ) _P .$$\n",
    "\n",
    "The change in enthalpy will be calculated with departure functions\n",
    "$$ \\Delta \\underline{H} = \\Delta \\underline{H}^\\mathrm{IG} + \\Delta [ \\underline{H} - \\underline{H}^\\mathrm{IG}] $$\n",
    "\n",
    "\n",
    "or\n",
    "\n",
    "$$\\Delta \\underbar{H} = \\Delta \\underline{H}^\\mathrm{IG} + [ \\underline{H} - \\underline{H}^\\mathrm{IG}]_2 - [ \\underline{H} - \\underline{H}^\\mathrm{IG}]_1.$$\n",
    "\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": 62,
   "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\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74031b60-3999-49f8-91a2-eebfeeb09d83",
   "metadata": {},
   "source": [
    "## Data for ethane\n",
    "\n",
    "For ethane, we need the critical pressure and temperature and Pitzer's acentric factor:\n",
    "\n",
    "$P_c = 4.884$ MPa   \n",
    "$T_c = 305.4$ K  \n",
    "$\\underline{V}_c = 0.184\\,\\mathrm{m^3/kmol} = 9.9\\times10^{-5}\\,\\mathrm{m^3/mol}$  \n",
    "$\\omega = 0.098$  \n",
    "\n",
    "We also need heat capacity data in the ideal state. We use the empirical equation in SIS A.II:\n",
    "\n",
    "$$ C^*_p(T) = A + BT + CT^2 + DT^3 $$\n",
    "\n",
    "The code block below has values for several molecules, including ethane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "7579e5f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Critical values for Ethane and the acentric factor\n",
    "# pressure in Pa, temperature in K\n",
    "Pc = 4.884e6\n",
    "Tc = 305.4\n",
    "omega = 0.098\n",
    "\n",
    "# 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[4,:]  # equation parameters for ethane"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79739b7d",
   "metadata": {},
   "source": [
    "## Function to calculate the heat capacity $C_p^*(T)$\n",
    "\n",
    "Here, the ideal state constant pressure heat capacity $ C_p^*(T)$ is given by an empirical equation\n",
    "\n",
    "$$ C_p^*(T) = a + bT + cT^2 + dT^3 \\tag{Table A.II}$$\n",
    "\n",
    "where $T$ is in absolute temperature Kelvins and $C_p^*$ is in units J/(mol K)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "7f48f67a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the heat capacity using SIS correlation in Appendix A.II\n",
    "\n",
    "def calc_Cp_IG(T, params):\n",
    "    A, B, C, D = params[:]\n",
    "    B *= 10**-2\n",
    "    C *= 10**-5\n",
    "    D *= 10**-9\n",
    "    return A + B*T + C*T**2 + D*T**3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b0c2408-14d8-4e9b-b593-414804278b6d",
   "metadata": {},
   "source": [
    "## Function to calculate the enthalpy change 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 SIS 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": 65,
   "id": "b587ff94-888c-455a-b216-f68a3a7ff462",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def calc_deltaH_IG(T1,T2,params):\n",
    "    A, B, C, D = params[:4]\n",
    "    B *= 10**-2\n",
    "    C *= 10**-5\n",
    "    D *= 10**-9\n",
    "    deltaHIG = A*(T2-T1) + B*(T2**2-T1**2)/2 + C*(T2**3-T1**3)/3 + D*(T2**4-T1**4)/4\n",
    "    return deltaHIG"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4133d8cd",
   "metadata": {},
   "source": [
    "## Function for Generalized Peng-Robinson departure function \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": 73,
   "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-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/(T*Tc)**0.5\n",
    "\n",
    "    depH = R*T*(Z-1)+(T*dadT-a)/(2*np.sqrt(2)*b) \\\n",
    "        *np.log((Z+(1+np.sqrt(2))*B)/(Z+(1-np.sqrt(2))*B))\n",
    "    return depH"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad238854-9629-4ce9-9f2a-0c92edcecc92",
   "metadata": {},
   "source": [
    "# Enthalpy calculation function\n",
    "\n",
    "We will define a function to calculate the enthalpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "327f42c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def del_H(T1, T2, P1, P2, params):\n",
    "    '''\n",
    "    params should include: A, B, C, D for Cp correlation    '''\n",
    "\n",
    "    # Calculate the compressibility factor\n",
    "    # Use the maximum value in case there are multiple roots\n",
    "    Z1 = np.max(PR_compressibility(P1,T1,Pc,Tc,omega))\n",
    "\n",
    "    # Calculate the first departure function\n",
    "    depH1 = calc_depH(T1,P1,Z1)\n",
    "\n",
    "    # Calculate the compressibility factor\n",
    "    # Use the maximum value in case there are multiple roots\n",
    "    Z2 = np.max(PR_compressibility(P2,T2,Pc,Tc,omega))\n",
    "\n",
    "    # Calculate the second departure function\n",
    "    depH2 = calc_depH(T2,P2,Z2)\n",
    "\n",
    "    # Calculate the ideal gas contribution\n",
    "    deltaH_IG = calc_deltaH_IG(T1,T2,params)\n",
    "\n",
    "    # Add for the total change in enthalpy\n",
    "    return deltaH_IG + (depH2 - depH1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dd02464",
   "metadata": {},
   "source": [
    "# Routine to calculate heat capacity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "bb8d3194-40ce-463c-9914-5eed37d4f7de",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Heat capacity for ethane at different pressures')"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 45360 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Known values in problem\n",
    "# The initial and final temperatures are\n",
    "Ti = 30 + 273.15     # Beginning temperature in K\n",
    "Tf = 200 + 273.15    # final temperature in K\n",
    "delta_T = 1 # use increments to calculate numerical derivative\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "for P in [0.1e6, 2e6, 4e6]:\n",
    "\n",
    "    T = np.arange(Ti, Tf + 1, delta_T)\n",
    "\n",
    "    cp_values = np.zeros_like(T, dtype=float)\n",
    "\n",
    "    for idx, T0 in enumerate(T):\n",
    "        # Compute ΔH over 1 K interval for each temperature\n",
    "        dH = del_H(T0, T0 + delta_T, P, P, i[:4])\n",
    "        cp_values[idx] = dH / delta_T\n",
    "\n",
    "    ax.plot(T, cp_values, label=f'P = {P/1e6:.1f} MPa')\n",
    "\n",
    "cp_ideal = calc_Cp_IG(T, i[:4])\n",
    "ax.plot(T, cp_ideal, '--', label='Ideal Gas')\n",
    "\n",
    "ax.legend()\n",
    "plt.xlabel('T (K)')\n",
    "plt.ylabel(\"$C_p$ (J/molK)\")\n",
    "plt.title('Heat capacity for ethane at different pressures')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5e87754",
   "metadata": {},
   "source": [
    "This makes sense. The isotherms on the pressure-enthalpy plot tend to bend toward lower values of enthalpy with higher pressures. That is, $(\\partial P/\\partial \\underbar{H})_T$ has a lower and negative slope on the plot. Thus, at $P = 4$ MPa, for each 10 degree difference in temperature, there is nearly double the amount of enthalpy change. Pretty interesting!\n",
    "\n",
    "Second, the non-lienar behavior also makes sense. This isobar eventually intersects the VLE curve (it’s below the critical pressure) and so you expect the heat capacity to diverge! Why? Think of the intersecting isotherm and what the slope represented by looks like in the VLE region. Then, as you go to higher temperatures along the isobar away from the VLE curve, you begin to recover more ideal-like behavior (and its temperature dependence)."
   ]
  }
 ],
 "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
}
