As in Michelsen's book, a generalized volume-translated cubic EOS can be given the common form: $$ p = \frac{RT}{v-b}-\frac{a(T)}{(v+\Delta_1b)(v+\Delta_2b)} $$
$$ p = \frac{RT}{v+c-b}-\frac{a(T)}{(v+c+\Delta_1b)(v+c+\Delta_2b)} $$
In [212]:
from __future__ import division
from sympy import *
from IPython.display import display, Math, Latex
from IPython.core.display import display_html
init_session(quiet=True)
init_printing()
In [213]:
rho_r,rho,tau,delta,T_r = symbols('rho_r,rho,tau,delta,T_r', positive = True, real = True)
R,Delta_1,Delta_2,T,v,c = symbols('R,Delta_1,Delta_2,T,v,c', positive = True, real = True)
x_i,x_j,x_k = symbols('x_i,x_j,x_k', positive=True, real = True)
a = symbols('a', cls=Function, positive = True, real = True)(tau)
b = symbols('b', positive = True, real = True)
In [214]:
normal_cubic = True # True for no volume translation
cubic = 'SRK'
In [215]:
p = R*T/(v+c-b)-a/((v+c+Delta_1*b)*(v+c+Delta_2*b))
if cubic == 'SRK':
p = p.subs(Delta_1,1).subs(Delta_2,0)
elif cubic == 'VdW':
p = p.subs(Delta_1,0).subs(Delta_2,0)
elif cubic == 'PR':
p = p.subs(Delta_1,1+sqrt(2)).subs(Delta_2,1-sqrt(2))
Math('p = '+latex(p))
Out[215]:
In [216]:
Z = p*v/(R*T)
Math('Z = '+latex(Z))
Out[216]:
In [217]:
Z = Z.replace(v, 1/rho)
display(Math('Z = '+latex(Z)))
Z = Z.replace(T, T_r/tau)
Z = Z.replace(rho, delta*rho_r)
display(Math('Z = '+latex(Z)))
dalphar_dDelta = (Z-1)/delta
dalphar_dDelta = expand(simplify(dalphar_dDelta))
Math('\\frac{d\\alpha^r}{d\delta}|\\tau = ' + latex(dalphar_dDelta))
Out[217]:
In [218]:
antideriv_pieces = 0
for arg in expand(simplify(dalphar_dDelta))._args:
antideriv_pieces += integrate(arg,delta)
antideriv_pieces = simplify(antideriv_pieces)
alphar = antideriv_pieces - antideriv_pieces.subs(delta,0)
display(Math('\\alpha^r = ' + latex(simplify(alphar))))
display(Math('\\lim_{c \\to 0}\\alpha^r = ' + latex(simplify(limit(alphar,c,0)))))
In [219]:
# Maybe turn off volume translation
if normal_cubic:
alphar = limit(alphar,c,0)
In [220]:
def format_deriv(arg, itau, idel, RHS):
"""
A function for giving a nice latex representation of
the partial derivative in question
"""
if itau+idel == 1:
numexp = ''
else:
numexp = '^{{{s:d}}}'.format(s=itau+idel)
if itau == 0:
tau = ''
elif itau == 1:
tau = '\\partial \\tau'
else:
tau = '\\partial \\tau^{{{s:d}}}'.format(s=itau)
if idel == 0:
delta = ''
elif idel == 1:
delta = '\\partial \\delta'
else:
delta = '\\partial \\delta^{{{s:d}}}'.format(s=idel)
temp = '\\frac{{\\partial{{{numexp:s}}} {arg:s}}}{{{{{tau:s}}}{{{delta:s}}}}} = '
if itau == 0 and idel == 0:
as_latex = arg +' = ' + latex(RHS)
else:
as_latex = temp.format(**locals()) + latex(RHS)
return Math(as_latex), as_latex
In [221]:
for deriv_count in range(1,5):
for dtau in range(deriv_count+1):
ddelta = deriv_count-dtau
#print dtau, ddelta
mth, tex = format_deriv('\\alpha^r', dtau, ddelta,
diff(diff(alphar,tau,dtau),delta,ddelta))
#print tex
display(mth)
In [222]:
a0,omega,kappa = symbols('a0,omega,kappa')
ai=a0*(1+kappa*(1-sqrt(1/tau)))**2
for diff_count in range(0,5):
mth, tex = format_deriv('a_i',diff_count, 0, diff(ai, tau, diff_count))
display(mth)