The iPSRV from Teus van der Stelt et al. $$\kappa = \kappa_0+\kappa_1\left\lbrace \sqrt{[A-D(T_r+B)]^2+E}+A-D(T_r+B) \right\rbrace\sqrt{T_r+C}$$ $$\kappa_0 = 0.378893 + 1.4897153\omega-0.17131848\omega^2+0.0196554\omega^3$$ $$p = \frac{RT}{v-b}-\frac{a}{v^2+2 b v-b^2}$$ with A = 1.1 B = 0.25 C = 0.2 D = 1.2 E = 0.01

For ethanol, the value $\kappa_1$ is -0.03374

For water, the value $\kappa_1$ is -0.06635



In [64]:

from CoolProp.CoolProp import Props
from math import sqrt
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
%matplotlib inline




In [65]:

kappa_dict = dict(Ethanol = -0.03374,
Water = -0.06635,
Propane = 0.03161)



## Pure Fluids



In [66]:

T = 300 #[K]
p = 3000 #[kPa]
Fluid = 'Propane'
rhoc = Props(Fluid,'rhocrit') #[kg/m^3]
omega = Props(Fluid,'accentric') #[-]
R = 8.314472/(Props(Fluid,'molemass'))
Tc = Props(Fluid,'Tcrit')
pc = Props(Fluid,'pcrit')
T_r = T/Tc
kappa_1 = kappa_dict[Fluid]
rho_EOS = Props('D','T',T,'P',p,Fluid)




In [67]:

Ts = np.linspace(Props(Fluid,'Ttriple'),Tc,300)
ps = Props('P','T',Ts,'Q',1,Fluid)
plt.plot(Ts,ps)
plt.xlabel('T [K]')
plt.ylabel('p [kPa]')
plt.yscale('log')
plt.plot(T,p,'o')




Out[67]:

[<matplotlib.lines.Line2D at 0x871c370>]




In [68]:

A = 1.1
B = 0.25
C = 0.2
D = 1.2
E = 0.01
kappa_0 = 0.378893 + 1.4897153*omega-0.17131848*omega**2+0.0196554*omega**3
kappa = kappa_0+kappa_1*(sqrt((A-D*(T_r+B))**2+E)+A-D*(T_r+B) )*sqrt(T_r+C)
alpha = (1+kappa*(1-sqrt(T_r)))**2
a = 0.457235*R**2*Tc**2/pc*alpha
b = 0.077796*R*Tc/pc




In [69]:

A = a*p/(R*R*T*T)
B = b*p/(R*T)




In [70]:

Z = np.linspace(0,1)
def f(Z):
return 1*Z**3+(-1+B)*Z**2+ (A-3*B*B-2*B)*Z -A*B+B**2+B**3
plt.plot(Z,f(Z))
plt.axhline(0,linestyle = '--')




Out[70]:

<matplotlib.lines.Line2D at 0x8909b50>




In [71]:

Z = np.roots([1,B-1,A-2*B-3*B*B,-A*B+B**2+B**3])
rho = p/(Z*R*T) # Z = p/(rho*R*T)
rho




Out[71]:

array([  69.42536782-63.64922759j,   69.42536782+63.64922759j,
518.88328530 +0.j        ])




In [72]:

rho_EOS




Out[72]:

495.5913139092507



# Mixtures



In [117]:

components = ['Ethanol','Water']
z = [0.4, 0.6]
#T = 562.8969961 #[K]
T = 400 #[K]
p = 12e6 #[Pa]




In [118]:

omega = np.array([Props(f,'accentric') for f in components]) #[-]
R = 8.314472 #[J/mol/K]
Tc = [Props(f,'Tcrit') for f in components] #[K]
pc = [Props(f,'pcrit')*1000 for f in components] #[Pa]

A = 1.1
B = 0.25
C = 0.2
D = 1.2
E = 0.01

kappa_0 = 0.378893 + 1.4897153*omega-0.17131848*omega**2+0.0196554*omega**3
kappa_1 = [kappa_dict[f] for f in components]

a,b = 0,0
for i in range(len(components)):

b += z[i]*0.077796*R*Tc[i]/pc[i]

for j in range(len(components)):
kappa_i = kappa_0[i]+kappa_1[i]*(sqrt((A-D*(T/Tc[i]+B))**2+E)+A-D*(T/Tc[i]+B) )*sqrt(T/Tc[i]+C)
alpha_i = (1+kappa_i*(1-sqrt(T/Tc[i])))**2
a_i = 0.457235*R**2*Tc[i]**2/pc[i]*alpha_i

kappa_j = kappa_0[j]+kappa_1[j]*(sqrt((A-D*(T/Tc[j]+B))**2+E)+A-D*(T/Tc[j]+B) )*sqrt(T/Tc[j]+C)
alpha_j = (1+kappa_j*(1-sqrt(T/Tc[j])))**2
a_j = 0.457235*R**2*Tc[j]**2/pc[j]*alpha_j

if i != j:
k_ij = -0.3
else:
k_ij = 0

a += z[i]*z[j]*(a_i*a_j)**0.5*(1-k_ij)

AA = a*p/(R**2*T**2)
BB = b*p/(R*T)




In [119]:

coeffs = [1,BB-1,AA-2*BB-3*BB**2,-AA*BB+BB**2+BB**3]




In [120]:

Z = np.linspace(0,1)
def f(Z):
return sum([a*Z**b for (a,b) in zip(coeffs,reversed(range(len(coeffs))))])
plt.plot(Z,f(Z))
plt.axhline(0,linestyle = '--')




Out[120]:

<matplotlib.lines.Line2D at 0x97cd2d0>




In [121]:

Z = np.roots(coeffs)
rho = p/(Z*R*T) # Z = p/(rho*R*T), rho in [mol/m^3]
rho




Out[121]:

array([  1234.16865590-3249.94868206j,   1234.16865590+3249.94868206j,
24827.37270034   +0.j        ])




In [ ]: