Last active
February 19, 2020 12:08
-
-
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
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"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