thermohub/chapter_6/PengRobinson_Me_throttle_example.ipynb
2026-03-16 08:32:38 -06:00

576 lines
18 KiB
Text
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

{
"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
}