{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7af16f2a",
   "metadata": {},
   "source": [
    "# Phase diagram calculations for carbon dioxide"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32d5dd49",
   "metadata": {},
   "source": [
    "Contruct a P-T phase diagram for carbon dioxide."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb1887c3",
   "metadata": {},
   "source": [
    "## Antoine equation\n",
    "\n",
    "The Antoine equation is an emprical equation of the form\n",
    "\n",
    "$$ \\log P^\\mathrm{sub}(T) = A - \\frac{B}{T + C} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "83c296bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Antoine equation given parameters A, B, and C\n",
    "# Pressure is returned in Pa\n",
    "def P_antoine(A, B, C, T):\n",
    "    return 10**(A - B/(T + C))*101325"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed693003",
   "metadata": {},
   "source": [
    "## Clausius-Clapeyron equation\n",
    "\n",
    "Given $\\Delta^\\mathrm{sat}\\underline{H}$, $T_1$, $P_1$, and $T_2$, return $P_2$ by the equation\n",
    "\n",
    "$$ \\ln \\frac{P^\\mathrm{sat}(T_2)}{P^\\mathrm{sat}(T_1)} = \\frac{-\\Delta^\\mathrm{sat}\\underline{H}}{R} \\left ( \\frac{1}{T_2} - \\frac{1}{T_1} \\right ) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4a95c76a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clausius-Clapeyron equation\n",
    "# Pressures in Pa, Temperatues in K\n",
    "# Enthalpy in J / mol\n",
    "\n",
    "import numpy as np\n",
    "from scipy import constants\n",
    "R = constants.R\n",
    "\n",
    "def P_CC(delta_sat_H, T1, P1, T2):\n",
    "    return P1*np.exp(-delta_sat_H/R*(1/T2-1/T1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93580680",
   "metadata": {},
   "source": [
    "## $P-T$ plot for CO2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87fec321",
   "metadata": {},
   "source": [
    "### Data for CO2 -- sublimation, VLE, critical point, and triple point\n",
    "\n",
    "NIST Webbook data is included below:\n",
    "\n",
    "https://webbook.nist.gov/cgi/cbook.cgi?ID=C124389&Units=SI&Mask=4#Thermo-Phase\n",
    "\n",
    "Full vapor-liquid properties can also be calculated in the NIST Webbook:\n",
    "\n",
    "https://webbook.nist.gov/cgi/fluid.cgi?TLow=217&THigh=304&TInc=&Digits=5&ID=C124389&Action=Load&Type=SatP&TUnit=K&PUnit=MPa&DUnit=mol%2Fm3&HUnit=kJ%2Fmol&WUnit=m%2Fs&VisUnit=uPa*s&STUnit=N%2Fm&RefState=DEF#Auxiliary\n",
    "\n",
    "We can find other thermodynamic data in:\n",
    "\n",
    "1. A. Jäger and R. Span, Equation of State for Solid Carbon Dioxide Based on the Gibbs Free Energy, J. Chem. Eng. Data 57, 590 (2012). https://doi.org/10.1021/je2011677\n",
    "1. R. Span and W. Wagner, A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa, Journal of Physical and Chemical Reference Data 25, 1509 (1996). https://doi.org/10.1063/1.555991\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2b4bdcfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data for CO2\n",
    "# Molecular weight\n",
    "MW = 44.0095  # g/mol\n",
    "\n",
    "# Critical parameters and acentricity\n",
    "Pc = 73.8*101325 # Critical pressure in Pa\n",
    "Tc = 304.18 # Critical temp in K\n",
    "Vc = 0.0000919 # Critical volume in m^3/mol\n",
    "omega = 0.225 # acentric factor\n",
    "\n",
    "# Triple point\n",
    "T_triple = 216.58 # temp in K\n",
    "P_triple = 5.185*101325 # pressure in Pa\n",
    "\n",
    "# Antoine equation parameters\n",
    "# Valid temperature range in K\n",
    "T_Antoine_min = 154.26 \n",
    "T_Antoine_max = 195.89\n",
    "\n",
    "A = 6.81228\n",
    "B = 1301.679\n",
    "C = -3.494\n",
    "\n",
    "# Enthalpy of sublimation \n",
    "H_sub_1 = 26.1*1000 # J / mol, 198 to 216 K\n",
    "H_sub_2 = 25.2*1000 # J / mol, 154 to 196 K\n",
    "\n",
    "# Enthalpy of vaporization\n",
    "H_vap_1 = 16.7*1000 # J / mol, 273 to 304 K\n",
    "H_vap_2 = 16.4*1000 # J / mol, 216 to 273 K\n",
    "H_vap_3 = 16.7*1000 # J / mol, 267 to 303 K"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72d6aac4",
   "metadata": {},
   "source": [
    "### Solid-liquid equilibrium\n",
    "\n",
    "Finding data for the solid-liquid equilibrium line is a little more challenging. Data for CO2 lists \n",
    "\n",
    "1562 kg/m3 solid at 1 atm (100 kPa) and −78.5 C \\\n",
    "1101 kg/m3 liquid at saturation −37 C \\\n",
    "1.977 kg/m3 gas at 1 atm (100 kPa) and 0 C\n",
    "\n",
    "Assuming that the solid and liquid specific (molar) volumes don't change much, we could estimate the fusion (melting) line with the Classius equation. But can we estimate the heat of fusion?\n",
    "\n",
    "J\\\"ager and Span (J. Chemical Eng. Data, 2006) cite the following value:\n",
    "\n",
    "$$\\Delta^\\mathrm{fus}\\underline{H} = 8875 \\text{ J/mol}$$\n",
    "\n",
    "By the Clapeyron equation, we can then find\n",
    "\n",
    "$$ \\left ( \\frac{\\partial P^\\mathrm{sat}}{\\partial T} \\right )_{\\underline{G}^I = \\underline{G}^{II}} = \\frac{\\Delta \\underline{H}}{T \\Delta \\underline{V} }$$\n",
    "\n",
    "As expected, this line should be quite steep. We can integrate the function to give\n",
    "\n",
    "$$ P^\\mathrm{sat}(T_2) = P^\\mathrm{sat}(T_1) + \\frac{\\Delta \\underline{H}}{\\Delta \\underline{V} } \\ln \\left ( \\frac{T_2}{T_1} \\right )$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "cd670958",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Integrated Clapeyron equation\n",
    "def P_clapeyron(T1, P1, T2, delta_H, delta_V):\n",
    "    return P1 + delta_H/delta_V*np.log(T2/T1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1f0615ef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.546529284164859e-05\n",
      "92965723.30973995\n"
     ]
    }
   ],
   "source": [
    "# Calculate Delta V\n",
    "delta_fus_V = 1/((1562-1101)/(MW/1000)) # Change in molar volume\n",
    "\n",
    "# Other data\n",
    "delta_fus_H = 8875 # Heat of fusion, J/mol\n",
    "\n",
    "print(delta_fus_V)\n",
    "print(delta_fus_H/delta_fus_V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3289cdae",
   "metadata": {},
   "source": [
    "### Building the PT diagram plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "69d95f54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'P-T diagram for carbon dioxide')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 43464 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Make a blank plot\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Plot the triple point\n",
    "plt.plot(T_triple, P_triple, \"bo\", label=\"triple point\")\n",
    "\n",
    "# Plot the critical point\n",
    "plt.plot(Tc, Pc, \"ro\", label=\"critical point\")\n",
    "\n",
    "# Plot the Antoine equation\n",
    "T = np.arange(T_Antoine_min, T_Antoine_max, 1)\n",
    "ax.plot(T, P_antoine(A, B, C, T), label=r\"$P^\\mathrm{sub}$ (Antoine)\")\n",
    "\n",
    "# Use Clausius-Clapeyron to plot the sublimation curve\n",
    "T = np.linspace(T_Antoine_max, T_triple, 300)\n",
    "P = P_CC(H_sub_1,T_triple, P_triple, T)\n",
    "ax.plot(T, P, label=r\"$P^\\mathrm{sub}$ (Clausius-Clapeyron)\")\n",
    "\n",
    "# Use Clausius-Clapeyron to plot the vaporization curve\n",
    "T = np.linspace(T_triple, Tc, 300)\n",
    "P = P_CC(H_vap_1,T_triple, P_triple, T)\n",
    "ax.plot(T, P, label=r\"$P^\\mathrm{vap}$ (Clausius-Clapeyron)\")\n",
    "\n",
    "# Use the integrated Clapeyron equation to plot the melting line\n",
    "T = np.linspace(T_triple, T_triple+18, 100)\n",
    "P = P_clapeyron(T_triple, P_triple, T, delta_fus_H, delta_fus_V)\n",
    "ax.plot(T, P, label=r\"$P^\\mathrm{fus}$ (Clausius)\")\n",
    "\n",
    "# Label and axes\n",
    "#plt.ylim([0,1e6])\n",
    "#ax.set_yscale(\"log\")\n",
    "ax.legend(frameon=True, fancybox=True, edgecolor='none')\n",
    "plt.xlabel('T (K)')\n",
    "plt.ylabel('P (Pa)')\n",
    "plt.title('P-T diagram for carbon dioxide')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa885660",
   "metadata": {},
   "source": [
    "This seems reasonable. You can see that the vaporization curve calculated by the Clausius-Clapeyron equation is senstive to the heat of vaporization. The value we use overestimates the critical point slightly. A slightly lower value would bring the two into agreement.\n",
    "\n",
    "Let's rebuild the plot, but focus on the lower part of the phase diagram by plotting it on a log-linear scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "0abc0fe5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'P-T diagram for carbon dioxide')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 42924 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Make a blank plot\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Plot the triple point\n",
    "plt.plot(T_triple, P_triple, \"bo\", label=\"triple point\")\n",
    "\n",
    "# Plot the critical point\n",
    "plt.plot(Tc, Pc, \"ro\", label=\"critical point\")\n",
    "\n",
    "# Plot the Antoine equation\n",
    "T = np.arange(T_Antoine_min, T_Antoine_max, 1)\n",
    "ax.plot(T, P_antoine(A, B, C, T), label=r\"$P^\\mathrm{sub}$ (Antoine)\")\n",
    "\n",
    "# Use Clausius-Clapeyron to plot the sublimation curve\n",
    "T = np.linspace(T_Antoine_max, T_triple, 300)\n",
    "P = P_CC(H_sub_1,T_triple, P_triple, T)\n",
    "ax.plot(T, P, label=r\"$P^\\mathrm{sub}$ (Clausius-Clapeyron)\")\n",
    "\n",
    "# Use Clausius-Clapeyron to plot the vaporization curve\n",
    "T = np.linspace(T_triple, Tc, 300)\n",
    "P = P_CC(H_vap_1,T_triple, P_triple, T)\n",
    "ax.plot(T, P, label=r\"$P^\\mathrm{vap}$ (Clausius-Clapeyron)\")\n",
    "\n",
    "# Use the integrated Clapeyron equation to plot the melting line\n",
    "T = np.linspace(T_triple, T_triple+18, 100)\n",
    "P = P_clapeyron(T_triple, P_triple, T, delta_fus_H, delta_fus_V)\n",
    "ax.plot(T, P, label=r\"$P^\\mathrm{fus}$ (Clausius)\")\n",
    "\n",
    "# Label and axes\n",
    "ax.set_yscale(\"log\")\n",
    "#plt.ylim([0,4e6])\n",
    "#plt.xlim([T_Antoine_min,240])\n",
    "ax.legend(frameon=True, fancybox=True, edgecolor='none')\n",
    "plt.xlabel('T (K)')\n",
    "plt.ylabel('P (Pa)')\n",
    "plt.title('P-T diagram for carbon dioxide')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26ddb75",
   "metadata": {},
   "source": [
    "## Answer for part C of homework problem\n",
    "\n",
    "In our exerpience with carbon dioxide, it is often in its gaseous state and sometimes in its solid form. The reason is clear: at atmospheric pressure, we are well below the triple point, the lowest temperature where we would observe a liquid carbon dioxide phase. \n",
    "\n",
    "When the gas phase is cooled (for instance, by a Joule-Thompson expansion) we expect to form the solid phase. Likewise, when handling dry ice, we observe that it sublimates directly to the gas phase without melting into a liquid first."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce26d98",
   "metadata": {},
   "source": [
    "## Answer for part B\n",
    "\n",
    "Part B asks us what the state of a full gas cylinder for a home carbonation system is at room temperature (298 K) given that the cylinder volume is 0.64 L and initially contains 410 g of\n",
    "carbon dioxide.\n",
    "\n",
    "In this case, we have the temperature and molar volume:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "f5a1a461",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Molar volume: 6.87e-05 m^3/mol\n"
     ]
    }
   ],
   "source": [
    "T = 298\n",
    "V_molar = 0.64/1000/(410/MW)\n",
    "\n",
    "print(f\"Molar volume: {V_molar:.2e} m^3/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d34536",
   "metadata": {},
   "source": [
    "From the temperature, we know that we are below the critical temperature, and the molar volume is much less than critical volume. It's reasonable to assume that the carbon dioxide is in a liquid phase at this temperature and molar volume. Looking at the PV diagram we calculated in problem 1, this is close to the saturated liquid. \n",
    "\n",
    "The maximum pressure on the cylinder is 250 bar, which is 25 MPa, but this is likely the limit of the cylinder pressure limit. The CO2 would be supercritical in this case.\n",
    "\n",
    "We can calculate the pressure using the Peng-Robinson equation of state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c3e0db15",
   "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",
    "import numpy as np\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",
    "# Calculate the molar volume 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_volume(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",
    "    # Convert real values of Z to molar volume\n",
    "    V = real_roots*R*T/P\n",
    "\n",
    "    return V\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5e07a484",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The pressure is: 6.54e+06 Pa\n"
     ]
    }
   ],
   "source": [
    "print(f\"The pressure is: {PR_pressure(V_molar,298,Pc,Tc,omega):.2e} Pa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc8959a7",
   "metadata": {},
   "source": [
    "The calculated pressure is very close to the vaporization curve.\n",
    "\n",
    "The pressure from the Clausius-Clapeyron equation at this temperature is a little higher:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "eecd9901",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P^vap(298K) = 6.62e+06 Pa\n"
     ]
    }
   ],
   "source": [
    "print(f\"P^vap(298K) = {P_CC(H_vap_1,T_triple, P_triple, 298):.2e} Pa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6eba2f5f",
   "metadata": {},
   "source": [
    "But we noted that the Clausius-Clapeyron equation might overestimate the vapor pressure curve. A more accurate calculation is necessary."
   ]
  }
 ],
 "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
}
