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)
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]:
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]:
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]:
In [72]:
rho_EOS
Out[72]:
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]:
In [121]:
Z = np.roots(coeffs)
rho = p/(Z*R*T) # Z = p/(rho*R*T), rho in [mol/m^3]
rho
Out[121]:
In [ ]: