Skip to content

Instantly share code, notes, and snippets.

@santiago-salas-v
Last active February 19, 2020 12:08
Show Gist options
  • Save santiago-salas-v/02d50da5ec2c4080cff329c6fdcbf1b8 to your computer and use it in GitHub Desktop.
Save santiago-salas-v/02d50da5ec2c4080cff329c6fdcbf1b8 to your computer and use it in GitHub Desktop.
Mechanical critical point, analytical solution according to Watson H. A., Barton, P. I. 2017. Reliable Flash Calculations: Part 3. A Nonsmooth Approach to Density Extrapolation and Pseudoproperty Evaluation, Ind. Eng. Chem. Res. 2017, 56, 14832−14847
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CEOS mechanical critical point (mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generic CEOS\n",
"\n",
"$P=\\frac{R T}{v-b}-\\frac{a}{(v+\\epsilon b)(v+\\sigma b)}$\n",
"\n",
"bzw.\n",
"\n",
"$Z=\\frac{v}{v-b}-\\frac{a}{(v+\\epsilon b)(v+\\sigma b)}\\frac{v}{R T}$\n",
"\n",
"$a=a(T,\\textbf{z}), b=b(\\textbf{z})$\n",
"\n",
"Mechanical critical point condition: $\\left(\\frac{d p}{d\\rho}\\right)_{T,x}=\\left(\\frac{d^2 p}{d\\rho^2}\\right)_{T,x}=0$\n",
"\n",
"Approach: eliminate a, RT and determine $v_{mc}=1/\\rho_{mc}$ such that condition is fulfilled"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from sympy import symbols, simplify, expand, diff, N, sqrt, root, collect\n",
"from sympy.solvers import solve\n",
"a,b,rt,p,epsilon,sigma,v,rho,R=symbols('a,b,rt,p,epsilon,sigma,v,rho,R')\n",
"p=rt/(v-b)-a/((v+sigma*b)*(v+epsilon*b))\n",
"dpdv_t=-rt/(v-b)**2+a/((v+sigma*b)*(v+epsilon*b))*(1/(v+epsilon*b)+1/(v+sigma*b))\n",
"#dpdrho_t=dpdv_t*(-v**2)\n",
"dpdrho_t=dpdv_t.subs(v,1/rho)*(-1/rho**2)\n",
"d2pdv2_t=2*rt/(v-b)**3-2*a/((v+epsilon*b)*(v+sigma*b))**2-2*a/((v+epsilon*b)*(v+sigma*b))*(1/(v+epsilon*b)**2+1/(v+sigma*b)**2)\n",
"#d2pdrho2_t=d2pdv2_t*v**4+dpdv_t*2*v**3\n",
"d2pdrho2_t=d2pdv2_t.subs(v,1/rho)*1/rho**4+dpdv_t.subs(v,1/rho)*2*1/rho**3\n",
"s1=solve(dpdrho_t,a)[0]\n",
"s2=solve(d2pdrho2_t,a)[0]\n",
"#soln_1=simplify(solve(d2pdrho2_t.subs(a,s1),v))\n",
"#soln_1=simplify(solve(d2pdrho2_t.subs(a,s1),rho))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check consistency"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0\n",
"0\n"
]
}
],
"source": [
"print(expand(dpdrho_t.subs(v,1/rho))-expand(diff(p.subs(v,1/rho),rho)))\n",
"print(expand(d2pdrho2_t.subs(v,1/rho))-expand(diff(dpdrho_t.subs(v,1/rho),rho)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get denominator from expression: $\\left(\\frac{d p}{d\\rho}\\right)_{T,x}=0$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{rt \\left(- b \\rho \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) \\left(b \\epsilon \\rho + b \\rho \\sigma + 2\\right) + \\left(b \\rho - 1\\right) \\left(b \\epsilon \\rho \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) + b \\rho \\sigma \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) - \\left(b \\epsilon \\rho + 1\\right)^{2} + \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) - \\left(b \\rho \\sigma + 1\\right)^{2}\\right)\\right)}{\\left(b \\rho - 1\\right)^{3} \\left(b \\epsilon \\rho \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) + b \\rho \\sigma \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) - \\left(b \\epsilon \\rho + 1\\right)^{2} + \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) - \\left(b \\rho \\sigma + 1\\right)^{2}\\right)}$"
],
"text/plain": [
"rt*(-b*rho*(b*epsilon*rho + 1)*(b*rho*sigma + 1)*(b*epsilon*rho + b*rho*sigma + 2) + (b*rho - 1)*(b*epsilon*rho*(b*epsilon*rho + 1)*(b*rho*sigma + 1) + b*rho*sigma*(b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*epsilon*rho + 1)**2 + (b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*rho*sigma + 1)**2))/((b*rho - 1)**3*(b*epsilon*rho*(b*epsilon*rho + 1)*(b*rho*sigma + 1) + b*rho*sigma*(b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*epsilon*rho + 1)**2 + (b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*rho*sigma + 1)**2))"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simplify(dpdrho_t.subs(a,s2))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"rt*(-b*rho*(b*epsilon*rho + 1)*(b*rho*sigma + 1)*(b*epsilon*rho + b*rho*sigma + 2) + (b*rho - 1)*(b*epsilon*rho*(b*epsilon*rho + 1)*(b*rho*sigma + 1) + b*rho*sigma*(b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*epsilon*rho + 1)**2 + (b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*rho*sigma + 1)**2))/((b*rho - 1)**3*(b*epsilon*rho*(b*epsilon*rho + 1)*(b*rho*sigma + 1) + b*rho*sigma*(b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*epsilon*rho + 1)**2 + (b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*rho*sigma + 1)**2))\n"
]
}
],
"source": [
"print(simplify(dpdrho_t.subs(a,s2)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get denominator from expression: $\\left(\\frac{d^2 p}{d\\rho^2}\\right)_{T,x}=0$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle - \\frac{2 rt \\left(\\left(b \\rho - 1\\right) \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) + \\left(b \\rho - 1\\right) \\left(\\left(b \\epsilon \\rho + 1\\right)^{2} + \\left(b \\rho \\sigma + 1\\right)^{2}\\right) + \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) \\left(b \\epsilon \\rho + b \\rho \\sigma + 2\\right)\\right)}{\\rho \\left(b \\rho - 1\\right)^{3} \\left(b \\epsilon \\rho + 1\\right) \\left(b \\rho \\sigma + 1\\right) \\left(b \\epsilon \\rho + b \\rho \\sigma + 2\\right)}$"
],
"text/plain": [
"-2*rt*((b*rho - 1)*(b*epsilon*rho + 1)*(b*rho*sigma + 1) + (b*rho - 1)*((b*epsilon*rho + 1)**2 + (b*rho*sigma + 1)**2) + (b*epsilon*rho + 1)*(b*rho*sigma + 1)*(b*epsilon*rho + b*rho*sigma + 2))/(rho*(b*rho - 1)**3*(b*epsilon*rho + 1)*(b*rho*sigma + 1)*(b*epsilon*rho + b*rho*sigma + 2))"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simplify(d2pdrho2_t.subs(a,s1))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-2*rt*((b*rho - 1)*(b*epsilon*rho + 1)*(b*rho*sigma + 1) + (b*rho - 1)*((b*epsilon*rho + 1)**2 + (b*rho*sigma + 1)**2) + (b*epsilon*rho + 1)*(b*rho*sigma + 1)*(b*epsilon*rho + b*rho*sigma + 2))/(rho*(b*rho - 1)**3*(b*epsilon*rho + 1)*(b*rho*sigma + 1)*(b*epsilon*rho + b*rho*sigma + 2))\n"
]
}
],
"source": [
"print(simplify(d2pdrho2_t.subs(a,s1)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## PR"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\rho_{mc}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Substitute EOS $\\sigma=1+\\sqrt{2}, \\epsilon=1-\\sqrt{2}$\n",
"* solve vor $v_{MC}=1/\\rho_{MC}$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"#expression_0 = -simplify(dpdrho_t.subs(a,s2))*(\n",
"# (b - v)**3*(b*epsilon*(b*epsilon + v)*(b*sigma + v) + b*sigma*(b*epsilon + v)*(b*sigma + v) - v*(b*epsilon + v)**2 + v*(b*epsilon + v)*(b*sigma + v) - v*(b*sigma + v)**2)\n",
"# )+simplify(d2pdrho2_t.subs(a,s1))*(\n",
"# (b - v)**3*(b*epsilon + v)*(b*sigma + v)*(b*epsilon + b*sigma + 2*v)\n",
"# )\n",
"expression_0 = -simplify(dpdrho_t.subs(a,s2))*(\n",
" (b*rho - 1)**3*(b*epsilon*rho*(b*epsilon*rho + 1)*(b*rho*sigma + 1) + b*rho*sigma*(b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*epsilon*rho + 1)**2 + (b*epsilon*rho + 1)*(b*rho*sigma + 1) - (b*rho*sigma + 1)**2)\n",
" )+simplify(d2pdrho2_t.subs(a,s1))*(\n",
" rho*(b*rho - 1)**3*(b*epsilon*rho + 1)*(b*rho*sigma + 1)*(b*epsilon*rho + b*rho*sigma + 2)\n",
" )\n",
"substitutions=[[sigma,1+sqrt(2)],[epsilon,1-sqrt(2)]]\n",
"#soln_2=solve((expression_0).subs(substitutions),v)\n",
"soln_2=solve(expression_0.subs(substitutions),rho)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{- \\frac{1}{3} - \\frac{2^{\\frac{2}{3}}}{3 \\sqrt[3]{4 + 3 \\sqrt{2}}} + \\frac{\\sqrt[3]{2} \\sqrt[3]{4 + 3 \\sqrt{2}}}{3}}{b}$"
],
"text/plain": [
"(-1/3 - 2**(2/3)/(3*(4 + 3*sqrt(2))**(1/3)) + 2**(1/3)*(4 + 3*sqrt(2))**(1/3)/3)/b"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#collect(simplify(1/soln_2[2]),3*b) # second root is real\n",
"rho_mc=collect(simplify(soln_2[2]),3*b) # second root is real\n",
"rho_mc"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(-1/3 - 2**(2/3)/(3*(4 + 3*sqrt(2))**(1/3)) + 2**(1/3)*(4 + 3*sqrt(2))**(1/3)/3)/b\n"
]
}
],
"source": [
"print(collect(simplify(soln_2[2]),3*b))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sample value with $b=8,4117\\cdot 10^{-5}~m^3/mol$\n",
"\n",
"Expected: $\\rho_{mc}=3.01~kmol/m^3$ (source: Watson H. A., Barton, P. I. 2017. Reliable Flash Calculations: Part 3. A Nonsmooth Approach to Density Extrapolation and Pseudoproperty Evaluation, Ind. Eng. Chem. Res. 2017, 56, 14832−14847)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 3008.62592034428$"
],
"text/plain": [
"3008.62592034428"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#N(1/soln_2[2].subs(b,8.4117e-5))\n",
"N(soln_2[2].subs(b,8.4117e-5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$T_{mc}$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"substitutions=[[rho,soln_2[2]],[sigma,1+sqrt(2)],[epsilon,1-sqrt(2)]]\n",
"rt_mc=solve(dpdrho_t.subs(substitutions),rt)[0]\n",
"t_mc=rt_mc/R"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{2 a \\left(- 12 \\cdot 2^{\\frac{5}{6}} \\sqrt[3]{4 + 3 \\sqrt{2}} - 16 \\sqrt[3]{8 + 6 \\sqrt{2}} - 15 \\sqrt{2} - 20 + 10 \\left(8 + 6 \\sqrt{2}\\right)^{\\frac{2}{3}} + 15 \\sqrt[6]{2} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}}\\right)}{R b \\left(- 39 \\cdot 2^{\\frac{5}{6}} \\sqrt[3]{4 + 3 \\sqrt{2}} - 52 \\sqrt[3]{2} \\sqrt[3]{4 + 3 \\sqrt{2}} + 48 \\sqrt{2} + 70 + 24 \\sqrt[6]{2} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}} + 19 \\cdot 2^{\\frac{2}{3}} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}}\\right)}$"
],
"text/plain": [
"2*a*(-12*2**(5/6)*(4 + 3*sqrt(2))**(1/3) - 16*(8 + 6*sqrt(2))**(1/3) - 15*sqrt(2) - 20 + 10*(8 + 6*sqrt(2))**(2/3) + 15*2**(1/6)*(4 + 3*sqrt(2))**(2/3))/(R*b*(-39*2**(5/6)*(4 + 3*sqrt(2))**(1/3) - 52*2**(1/3)*(4 + 3*sqrt(2))**(1/3) + 48*sqrt(2) + 70 + 24*2**(1/6)*(4 + 3*sqrt(2))**(2/3) + 19*2**(2/3)*(4 + 3*sqrt(2))**(2/3)))"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t_mc"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2*a*(-12*2**(5/6)*(4 + 3*sqrt(2))**(1/3) - 16*(8 + 6*sqrt(2))**(1/3) - 15*sqrt(2) - 20 + 10*(8 + 6*sqrt(2))**(2/3) + 15*2**(1/6)*(4 + 3*sqrt(2))**(2/3))/(R*b*(-39*2**(5/6)*(4 + 3*sqrt(2))**(1/3) - 52*2**(1/3)*(4 + 3*sqrt(2))**(1/3) + 48*sqrt(2) + 70 + 24*2**(1/6)*(4 + 3*sqrt(2))**(2/3) + 19*2**(2/3)*(4 + 3*sqrt(2))**(2/3)))\n"
]
}
],
"source": [
"print(t_mc)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{0.17014442007035 a}{R b}$"
],
"text/plain": [
"0.17014442007035*a/(R*b)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(t_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$p_{mc}$ , directly from EOS at $T_{mc},\\rho_{mc}$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{a \\left(-760 - 537 \\sqrt{2} - 119 \\sqrt[3]{8 + 6 \\sqrt{2}} - 84 \\cdot 2^{\\frac{5}{6}} \\sqrt[3]{4 + 3 \\sqrt{2}} + 231 \\sqrt[6]{2} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}} + 164 \\left(8 + 6 \\sqrt{2}\\right)^{\\frac{2}{3}}\\right)}{2 b^{2} \\left(-1396 - 987 \\sqrt{2} + 69 \\sqrt[6]{2} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}} + 50 \\cdot 2^{\\frac{2}{3}} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}} + 276 \\cdot 2^{\\frac{5}{6}} \\sqrt[3]{4 + 3 \\sqrt{2}} + 391 \\sqrt[3]{2} \\sqrt[3]{4 + 3 \\sqrt{2}}\\right)}$"
],
"text/plain": [
"a*(-760 - 537*sqrt(2) - 119*(8 + 6*sqrt(2))**(1/3) - 84*2**(5/6)*(4 + 3*sqrt(2))**(1/3) + 231*2**(1/6)*(4 + 3*sqrt(2))**(2/3) + 164*(8 + 6*sqrt(2))**(2/3))/(2*b**2*(-1396 - 987*sqrt(2) + 69*2**(1/6)*(4 + 3*sqrt(2))**(2/3) + 50*2**(2/3)*(4 + 3*sqrt(2))**(2/3) + 276*2**(5/6)*(4 + 3*sqrt(2))**(1/3) + 391*2**(1/3)*(4 + 3*sqrt(2))**(1/3)))"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"substitutions=[[v,1/soln_2[2]],[rt,R*t_mc],[sigma,1+sqrt(2)],[epsilon,1-sqrt(2)]]\n",
"p_mc = p.subs(substitutions)\n",
"simplify(p_mc)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{0.0132365678781271 a}{b^{2}}$"
],
"text/plain": [
"0.0132365678781271*a/b**2"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(p_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$Z_{mc}$ from definition $Z=\\frac{p v}{R T}$"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{- 1908 \\sqrt[6]{2} - 1349 \\cdot 2^{\\frac{2}{3}} + 84 \\cdot 2^{\\frac{5}{6}} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}} + 119 \\sqrt[3]{2} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}} + 537 \\sqrt{2} \\sqrt[3]{4 + 3 \\sqrt{2}} + 760 \\sqrt[3]{4 + 3 \\sqrt{2}}}{4 \\left(- 1246 \\sqrt[6]{2} - 881 \\cdot 2^{\\frac{2}{3}} + 185 \\sqrt{2} \\sqrt[3]{4 + 3 \\sqrt{2}} + 262 \\sqrt[3]{4 + 3 \\sqrt{2}} + 120 \\cdot 2^{\\frac{5}{6}} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}} + 170 \\sqrt[3]{2} \\left(4 + 3 \\sqrt{2}\\right)^{\\frac{2}{3}}\\right)}$"
],
"text/plain": [
"(-1908*2**(1/6) - 1349*2**(2/3) + 84*2**(5/6)*(4 + 3*sqrt(2))**(2/3) + 119*2**(1/3)*(4 + 3*sqrt(2))**(2/3) + 537*sqrt(2)*(4 + 3*sqrt(2))**(1/3) + 760*(4 + 3*sqrt(2))**(1/3))/(4*(-1246*2**(1/6) - 881*2**(2/3) + 185*sqrt(2)*(4 + 3*sqrt(2))**(1/3) + 262*(4 + 3*sqrt(2))**(1/3) + 120*2**(5/6)*(4 + 3*sqrt(2))**(2/3) + 170*2**(1/3)*(4 + 3*sqrt(2))**(2/3)))"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z_mc = p_mc * 1/rho_mc / (R * t_mc)\n",
"simplify(z_mc)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 0.307401308698703$"
],
"text/plain": [
"0.307401308698703"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(z_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## SRK"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\rho_{mc}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Substitute EOS $\\sigma=1, \\epsilon=0$\n",
"* solve vor $v_{MC}=1/\\rho_{MC}$"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"substitutions=[[sigma,1],[epsilon,0]]\n",
"#soln_2=solve((expression_0).subs(substitutions),v)\n",
"soln_2=solve(expression_0.subs(substitutions),rho)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{-1 + \\sqrt[3]{2}}{b}$"
],
"text/plain": [
"(-1 + 2**(1/3))/b"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#collect(simplify(1/soln_2[2]),3*b) # second root is real\n",
"rho_mc=collect(simplify(soln_2[0]),3*b) # first root is real\n",
"rho_mc"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(-1 + 2**(1/3))/b\n"
]
}
],
"source": [
"print(collect(simplify(soln_2[0]),3*b))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sample value with $b=8,4117\\cdot 10^{-5}~m^3/mol$\n",
"\n",
"Expected: $\\rho_{mc}=3.01~kmol/m^3$ (source: Watson H. A., Barton, P. I. 2017. Reliable Flash Calculations: Part 3. A Nonsmooth Approach to Density Extrapolation and Pseudoproperty Evaluation, Ind. Eng. Chem. Res. 2017, 56, 14832−14847)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 3089.99429241263$"
],
"text/plain": [
"3089.99429241263"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#N(1/soln_2[2].subs(b,8.4117e-5))\n",
"N(soln_2[0].subs(b,8.4117e-5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$T_{mc}$"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"substitutions=[[rho,soln_2[0]],[sigma,1],[epsilon,0]]\n",
"rt_mc=solve(dpdrho_t.subs(substitutions),rt)[0]\n",
"t_mc=rt_mc/R"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{a \\left(2 - \\sqrt[3]{2}\\right)^{3}}{2 R b}$"
],
"text/plain": [
"a*(2 - 2**(1/3))**3/(2*R*b)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t_mc"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a*(2 - 2**(1/3))**3/(2*R*b)\n"
]
}
],
"source": [
"print(t_mc)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{0.202676856535359 a}{R b}$"
],
"text/plain": [
"0.202676856535359*a/(R*b)"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(t_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$p_{mc}$ , directly from EOS at $T_{mc},\\rho_{mc}$"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{a \\left(- 3 \\cdot 2^{\\frac{2}{3}} + 1 + 3 \\sqrt[3]{2}\\right)}{b^{2}}$"
],
"text/plain": [
"a*(-3*2**(2/3) + 1 + 3*2**(1/3))/b**2"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"substitutions=[[v,1/soln_2[0]],[rt,R*t_mc],[sigma,1],[epsilon,0]]\n",
"p_mc = p.subs(substitutions)\n",
"simplify(p_mc)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{0.0175599937800211 a}{b^{2}}$"
],
"text/plain": [
"0.0175599937800211*a/b**2"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(p_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$Z_{mc}$ from definition $Z=\\frac{p v}{R T}$"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{1}{3}$"
],
"text/plain": [
"1/3"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z_mc = p_mc * 1/rho_mc / (R * t_mc)\n",
"simplify(z_mc)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 0.333333333333333$"
],
"text/plain": [
"0.333333333333333"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(z_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## VdW"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\rho_{mc}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Substitute EOS $\\sigma=0, \\epsilon=0$\n",
"* solve vor $v_{MC}=1/\\rho_{MC}$"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"substitutions=[[sigma,0],[epsilon,0]]\n",
"#soln_2=solve((expression_0).subs(substitutions),v)\n",
"soln_2=solve(expression_0.subs(substitutions),rho)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{1}{3 b}$"
],
"text/plain": [
"1/(3*b)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#collect(simplify(1/soln_2[2]),3*b) # second root is real\n",
"rho_mc=collect(simplify(soln_2[0]),3*b) # first root is real\n",
"rho_mc"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1/(3*b)\n"
]
}
],
"source": [
"print(collect(simplify(soln_2[0]),3*b))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sample value with $b=8,4117\\cdot 10^{-5}~m^3/mol$\n",
"\n",
"Expected: $\\rho_{mc}=3.01~kmol/m^3$ (source: Watson H. A., Barton, P. I. 2017. Reliable Flash Calculations: Part 3. A Nonsmooth Approach to Density Extrapolation and Pseudoproperty Evaluation, Ind. Eng. Chem. Res. 2017, 56, 14832−14847)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 3962.73444527662$"
],
"text/plain": [
"3962.73444527662"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#N(1/soln_2[2].subs(b,8.4117e-5))\n",
"N(soln_2[0].subs(b,8.4117e-5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$T_{mc}$"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"substitutions=[[rho,soln_2[0]],[sigma,0],[epsilon,0]]\n",
"rt_mc=solve(dpdrho_t.subs(substitutions),rt)[0]\n",
"t_mc=rt_mc/R"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{8 a}{27 R b}$"
],
"text/plain": [
"8*a/(27*R*b)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t_mc"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8*a/(27*R*b)\n"
]
}
],
"source": [
"print(t_mc)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{0.296296296296296 a}{R b}$"
],
"text/plain": [
"0.296296296296296*a/(R*b)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(t_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$p_{mc}$ , directly from EOS at $T_{mc},\\rho_{mc}$"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{a}{27 b^{2}}$"
],
"text/plain": [
"a/(27*b**2)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"substitutions=[[v,1/soln_2[0]],[rt,R*t_mc],[sigma,0],[epsilon,0]]\n",
"p_mc = p.subs(substitutions)\n",
"simplify(p_mc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$Z_{mc}$ from definition $Z=\\frac{p v}{R T}$"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{3}{8}$"
],
"text/plain": [
"3/8"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z_mc = p_mc * 1/rho_mc / (R * t_mc)\n",
"simplify(z_mc)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 0.375$"
],
"text/plain": [
"0.375000000000000"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(z_mc)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment