517 lines
15 KiB
Text
517 lines
15 KiB
Text
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "f225380f",
|
||
"metadata": {},
|
||
"source": [
|
||
"# Homework: Calculate the temperature of ethane in a throttling process\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 ethane through a throttling process.\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": 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",
|
||
"Ethane gas undergoes a rapid, adiabatic expansion through a continous throttling process from upstream conditions ${\\text{100}^\\circ}\\text{C}$ and 100 bar to atmospheric pressure.\n",
|
||
"\n",
|
||
"Calculate the downstream gas temperature.\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 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$ "
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"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",
|
||
"# Known values in problem\n",
|
||
"T1 = 373 # Inlet temperature in K\n",
|
||
"P1 = 10.e6 # Inlet pressure in Pa\n",
|
||
"P2 = 101325 # Outlet temperature in Pa"
|
||
]
|
||
},
|
||
{
|
||
"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": "43ba6036",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"1.2213490504256714 2.0475020475020473\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"TR1 = T1/Tc # Reduced temperature\n",
|
||
"PR1 = P1/Pc # Reduced pressure\n",
|
||
"\n",
|
||
"print(TR1,PR1)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"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": 6,
|
||
"id": "ca49e458-2711-465b-add8-198d80ff5131",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Z1 = 0.61\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"print(f\"Z1 = {Z1:.2f}\")"
|
||
]
|
||
},
|
||
{
|
||
"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": 7,
|
||
"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": "markdown",
|
||
"id": "c8640217",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Calculate the inlet departure function $[\\underbar{H}-\\underbar{H}^\\mathrm{IG}]_1$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 8,
|
||
"id": "c3dc8ef7-f035-4079-8e90-bb4cc500fbf9",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"*** depH1 = -5108 J/mol\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"depH1 = calc_depH(T1,P1,Z1)\n",
|
||
"\n",
|
||
"print(f\"*** depH1 = {depH1:.0f} J/mol\")"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "1fb6ec34",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Guess an outlet temperature\n",
|
||
"\n",
|
||
"Now, the challenge here is that the outlet temperature is not known. So, we have to guess an outlet temperature and calculate both the ideal contributino and the outlet departure function value."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 9,
|
||
"id": "02096ffe",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# GUESS OUTLET TEMP\n",
|
||
"T2 = 284.4 # outlet temp in K\n",
|
||
"T2R = T2/Tc # reduced outlet temperature"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "9b0c2408-14d8-4e9b-b593-414804278b6d",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Function to calculate the enthalpy change in the ideal state $\\Delta \\underbar{H}^\\mathrm{IG}$\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": 10,
|
||
"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[4,:] # equation parameters for ethane\n",
|
||
"\n",
|
||
"deltaHIG = i[0]*(T2-T1) + i[1]*10**-2*(T2**2-T1**2)/2 + i[2]*10**-5*(T2**3-T1**3)/3 + i[3]*10**-9*(T2**4-T1**4)/4"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "2dd02464",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Calculate the enthalpy change of the ideal state"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 11,
|
||
"id": "bb8d3194-40ce-463c-9914-5eed37d4f7de",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"*** Delta HIG = -5043 J/mol\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"print(\"*** Delta HIG = {:.0f} J/mol\".format(deltaHIG))"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "a5cf19a3-a21e-4fc0-a034-30d48834da91",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Calculate the outlet departure function $[\\underbar{H}-\\underbar{H}^\\mathrm{IG}]_2$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 12,
|
||
"id": "04bda783-43d6-4fb6-8f69-e87dc1a5a2be",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"depH2 = -63 J/mol\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"# 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(f\"depH2 = {depH2:.0f} J/mol\")\n",
|
||
"# , deltaH = {:.0f} J/mol <-- should be zero\".format(depH2,deltaHIG+depH2-depH1))"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "b59d01f8",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Check if the sum is zero $\\Delta \\underbar{H}^\\mathrm{IG} + [\\underbar{H}-\\underbar{H}^\\mathrm{IG}]_2 - [\\underbar{H}-\\underbar{H}^\\mathrm{IG}]_1 = 0$\n",
|
||
"Since $\\Delta \\underbar{H} = 0$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 13,
|
||
"id": "b30ecfe0",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"delta_H = 2 J/mol\n",
|
||
"for T2 = 284.40 K\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"print(f\"delta_H = {deltaHIG + depH2 - depH1:.0f} J/mol\")\n",
|
||
"print(f\"for T2 = {T2:.2f} K\")"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "8044d9eb",
|
||
"metadata": {},
|
||
"source": [
|
||
"If this doesn't equal zero (or is close to zero), guess another $T_2$ and recalculate $\\Delta \\underbar{H}^\\mathrm{IG}$ and $[\\underbar{H} - \\underbar{H}^\\mathrm{IG}]_2$."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "edee9974",
|
||
"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
|
||
}
|