Helmholtz energy conversion of cubics

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()


IPython console for SymPy 0.7.6.1 (Python 2.7.10-32-bit) (ground types: python)

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]:
$$p = \frac{R T}{- b + c + v} - \frac{a{\left (\tau \right )}}{\left(c + v\right) \left(b + c + v\right)}$$

In [216]:
Z = p*v/(R*T)
Math('Z = '+latex(Z))


Out[216]:
$$Z = \frac{v}{R T} \left(\frac{R T}{- b + c + v} - \frac{a{\left (\tau \right )}}{\left(c + v\right) \left(b + c + v\right)}\right)$$

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))


$$Z = \frac{1}{R T \rho} \left(\frac{R T}{- b + c + \frac{1}{\rho}} - \frac{a{\left (\tau \right )}}{\left(c + \frac{1}{\rho}\right) \left(b + c + \frac{1}{\rho}\right)}\right)$$
$$Z = \frac{\tau}{R T_{r} \delta \rho_{r}} \left(\frac{R T_{r}}{\tau \left(- b + c + \frac{1}{\delta \rho_{r}}\right)} - \frac{a{\left (\tau \right )}}{\left(c + \frac{1}{\delta \rho_{r}}\right) \left(b + c + \frac{1}{\delta \rho_{r}}\right)}\right)$$
Out[217]:
$$\frac{d\alpha^r}{d\delta}|\tau = - \frac{\tau a{\left (\tau \right )}}{R T_{r} b c \delta^{2} \rho_{r} + R T_{r} b \delta + R T_{r} c^{2} \delta^{2} \rho_{r} + 2 R T_{r} c \delta + \frac{R T_{r}}{\rho_{r}}} + \frac{\tau}{- b \delta^{2} \rho_{r} \tau + c \delta^{2} \rho_{r} \tau + \delta \tau} - \frac{1}{\delta}$$

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)))))


$$\alpha^r = \frac{1}{R T_{r} b} \left(R T_{r} b \log{\left (- \frac{1}{\rho_{r} \left(b - c\right)} \right )} - R T_{r} b \log{\left (\frac{\delta \rho_{r} \left(b - c\right) - 1}{\rho_{r} \left(b - c\right)} \right )} + a{\left (\tau \right )} \log{\left (\frac{c^{\tau} \left(b + c \delta \rho_{r} \left(b + c\right) + c\right)^{\tau}}{\left(c \left(\delta \rho_{r} \left(b + c\right) + 1\right)\right)^{\tau} \left(b + c\right)^{\tau}} \right )}\right)$$
$$\lim_{c \to 0}\alpha^r = - \log{\left (b \delta \rho_{r} - 1 \right )} + i \pi - \frac{\tau a{\left (\tau \right )}}{R T_{r} b} \log{\left (b \delta \rho_{r} + 1 \right )}$$

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)


$$\frac{\partial{} \alpha^r}{{}{\partial \delta}} = \frac{1}{R T_{r} b} \left(- \frac{R T_{r} b}{\delta - \frac{1}{b \rho_{r}}} - \frac{b \rho_{r} \tau a{\left (\tau \right )}}{b \delta \rho_{r} + 1}\right)$$
$$\frac{\partial{} \alpha^r}{{\partial \tau}{}} = \frac{1}{R T_{r} b} \left(- \tau \log{\left (b \delta \rho_{r} + 1 \right )} \frac{d}{d \tau} a{\left (\tau \right )} - a{\left (\tau \right )} \log{\left (b \delta \rho_{r} + 1 \right )}\right)$$
$$\frac{\partial{^{2}} \alpha^r}{{}{\partial \delta^{2}}} = \frac{1}{R T_{r}} \left(\frac{R T_{r}}{\left(\delta - \frac{1}{b \rho_{r}}\right)^{2}} + \frac{b \rho_{r}^{2} \tau a{\left (\tau \right )}}{\left(b \delta \rho_{r} + 1\right)^{2}}\right)$$
$$\frac{\partial{^{2}} \alpha^r}{{\partial \tau}{\partial \delta}} = \frac{1}{R T_{r} b} \left(- \frac{b \rho_{r} \tau \frac{d}{d \tau} a{\left (\tau \right )}}{b \delta \rho_{r} + 1} - \frac{b \rho_{r} a{\left (\tau \right )}}{b \delta \rho_{r} + 1}\right)$$
$$\frac{\partial{^{2}} \alpha^r}{{\partial \tau^{2}}{}} = - \frac{1}{R T_{r} b} \left(\tau \frac{d^{2}}{d \tau^{2}} a{\left (\tau \right )} + 2 \frac{d}{d \tau} a{\left (\tau \right )}\right) \log{\left (b \delta \rho_{r} + 1 \right )}$$
$$\frac{\partial{^{3}} \alpha^r}{{}{\partial \delta^{3}}} = - \frac{1}{R T_{r}} \left(\frac{2 R T_{r}}{\left(\delta - \frac{1}{b \rho_{r}}\right)^{3}} + \frac{2 b^{2} \rho_{r}^{3} \tau a{\left (\tau \right )}}{\left(b \delta \rho_{r} + 1\right)^{3}}\right)$$
$$\frac{\partial{^{3}} \alpha^r}{{\partial \tau}{\partial \delta^{2}}} = \frac{b \rho_{r}^{2} \left(\tau \frac{d}{d \tau} a{\left (\tau \right )} + a{\left (\tau \right )}\right)}{R T_{r} \left(b \delta \rho_{r} + 1\right)^{2}}$$
$$\frac{\partial{^{3}} \alpha^r}{{\partial \tau^{2}}{\partial \delta}} = - \frac{\rho_{r}}{R T_{r} \left(b \delta \rho_{r} + 1\right)} \left(\tau \frac{d^{2}}{d \tau^{2}} a{\left (\tau \right )} + 2 \frac{d}{d \tau} a{\left (\tau \right )}\right)$$
$$\frac{\partial{^{3}} \alpha^r}{{\partial \tau^{3}}{}} = - \frac{1}{R T_{r} b} \left(\tau \frac{d^{3}}{d \tau^{3}} a{\left (\tau \right )} + 3 \frac{d^{2}}{d \tau^{2}} a{\left (\tau \right )}\right) \log{\left (b \delta \rho_{r} + 1 \right )}$$
$$\frac{\partial{^{4}} \alpha^r}{{}{\partial \delta^{4}}} = \frac{1}{R T_{r}} \left(\frac{6 R T_{r}}{\left(\delta - \frac{1}{b \rho_{r}}\right)^{4}} + \frac{6 b^{3} \rho_{r}^{4} \tau a{\left (\tau \right )}}{\left(b \delta \rho_{r} + 1\right)^{4}}\right)$$
$$\frac{\partial{^{4}} \alpha^r}{{\partial \tau}{\partial \delta^{3}}} = - \frac{2 b^{2} \rho_{r}^{3}}{R T_{r} \left(b \delta \rho_{r} + 1\right)^{3}} \left(\tau \frac{d}{d \tau} a{\left (\tau \right )} + a{\left (\tau \right )}\right)$$
$$\frac{\partial{^{4}} \alpha^r}{{\partial \tau^{2}}{\partial \delta^{2}}} = \frac{b \rho_{r}^{2}}{R T_{r} \left(b \delta \rho_{r} + 1\right)^{2}} \left(\tau \frac{d^{2}}{d \tau^{2}} a{\left (\tau \right )} + 2 \frac{d}{d \tau} a{\left (\tau \right )}\right)$$
$$\frac{\partial{^{4}} \alpha^r}{{\partial \tau^{3}}{\partial \delta}} = - \frac{\rho_{r}}{R T_{r} \left(b \delta \rho_{r} + 1\right)} \left(\tau \frac{d^{3}}{d \tau^{3}} a{\left (\tau \right )} + 3 \frac{d^{2}}{d \tau^{2}} a{\left (\tau \right )}\right)$$
$$\frac{\partial{^{4}} \alpha^r}{{\partial \tau^{4}}{}} = - \frac{1}{R T_{r} b} \left(\tau \frac{d^{4}}{d \tau^{4}} a{\left (\tau \right )} + 4 \frac{d^{3}}{d \tau^{3}} a{\left (\tau \right )}\right) \log{\left (b \delta \rho_{r} + 1 \right )}$$

Derivatives of $a_i(\tau)$


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)


$$a_i = a_{0} \left(\kappa \left(1 - \frac{1}{\sqrt{\tau}}\right) + 1\right)^{2}$$
$$\frac{\partial{} a_i}{{\partial \tau}{}} = \frac{a_{0} \kappa}{\tau^{\frac{3}{2}}} \left(\kappa \left(1 - \frac{1}{\sqrt{\tau}}\right) + 1\right)$$
$$\frac{\partial{^{2}} a_i}{{\partial \tau^{2}}{}} = \frac{a_{0} \kappa}{2} \left(\frac{\kappa}{\tau^{3}} - \frac{1}{\tau^{\frac{5}{2}}} \left(3 \kappa \left(1 - \frac{1}{\sqrt{\tau}}\right) + 3\right)\right)$$
$$\frac{\partial{^{3}} a_i}{{\partial \tau^{3}}{}} = \frac{3 a_{0}}{4} \kappa \left(- \frac{3 \kappa}{\tau^{4}} + \frac{1}{\tau^{\frac{7}{2}}} \left(5 \kappa \left(1 - \frac{1}{\sqrt{\tau}}\right) + 5\right)\right)$$
$$\frac{\partial{^{4}} a_i}{{\partial \tau^{4}}{}} = \frac{3 a_{0}}{8} \kappa \left(\frac{29 \kappa}{\tau^{5}} - \frac{1}{\tau^{\frac{9}{2}}} \left(35 \kappa \left(1 - \frac{1}{\sqrt{\tau}}\right) + 35\right)\right)$$