{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3c29e0e8",
   "metadata": {},
   "source": "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/url/https%3A%2F%2Flem.che.udel.edu%2Fgit%2Ffurst%2Fthermohub%2Fraw%2Fbranch%2Fmain%2Fmodule_6%2FPeng_Robinson_EOS.ipynb)\n\n# Generalized Peng-Robinson Equation of State\n\nRoutines to calcualte the Generalized Peng-Robinson Equation of State\n\nSIS is Stanley I. Sandler, *Chemical, Biochemcial and Engineering Thermodynamics*, 5th ed.\n\nEric Furst  \nNovember 2025\n\nThe 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\nwith\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\nThe acentric factor $\\omega$ and the crticial temperatures and pressures are given in SIS table 6.6-1.\n\nCalculating 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\nwhere $Z$ is the compressibility factor\n\n$$ Z = \\frac{P \\underbar{V{}}}{RT} $$\n\nFor the Peng-Robinson EOS (see SIS Table 6.4-3),\n\n$$ \\alpha = -1 + B $$\n$$ \\beta = A - 3B^2 -2B $$\n$$ \\gamma = -AB + B^2 + B^3 $$\n\nand \n\n$$ A = \\frac{aP}{(RT)^2} $$\n$$ B = \\frac{bP}{RT} $$"
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "47da809f",
   "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": "markdown",
   "id": "2f55b187",
   "metadata": {},
   "source": [
    "## Example calculation\n",
    "\n",
    "Given the PR EOS, calculate the molar volume(s) of ethane given different temperature and pressures for ethane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "fd028ebd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6.08349837e-05 2.38749240e-04 1.65668503e-03]\n",
      "[494.28796012 125.94804501  18.15070421]\n"
     ]
    }
   ],
   "source": [
    "# Given the PR EOS, calculate the molar volume(s) of ethane given different T and P:\n",
    "\n",
    "# Data for ethane\n",
    "Pc = 4.884e6\n",
    "Tc = 305.4\n",
    "omega = 0.098\n",
    "MW = 30.07\n",
    "\n",
    "P = 1e6\n",
    "T = -33+273.15\n",
    "\n",
    "# molar volume m^3/mol\n",
    "print((PR_volume(P, T, Pc, Tc, omega))) \n",
    "\n",
    "# specific density kg/m^3\n",
    "print(1/(PR_volume(P, T, Pc, Tc, omega)/MW*1000)) \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f387f98f",
   "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
}