mirror of
https://github.com/CoolProp/CoolProp.git
synced 2026-01-12 15:38:29 -05:00
6913fc41dc
``` find . -regextype posix-extended -regex '.*\.(cpp|hpp|c|h|cxx|hxx|py)$' | xargs -I@ sed -i 's/[ \t]*$//' "@" ```
190 lines
8.3 KiB
Python
190 lines
8.3 KiB
Python
from __future__ import print_function
|
|
|
|
from math import sqrt, exp
|
|
import CoolProp
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import scipy.optimize
|
|
from math import log
|
|
|
|
|
|
def viscosity_dilute(fluid, T, e_k, sigma):
|
|
"""
|
|
T in [K], e_K in [K], sigma in [nm]
|
|
viscosity returned is in [Pa-s]
|
|
"""
|
|
Tstar = T / e_k
|
|
molemass = CoolProp.CoolProp.PropsSI(fluid, 'molemass') * 1000
|
|
|
|
# From Neufeld, 1972, Journal of Chemical Physics - checked coefficients
|
|
OMEGA_2_2 = 1.16145 * pow(Tstar, -0.14874) + 0.52487 * exp(-0.77320 * Tstar) + 2.16178 * exp(-2.43787 * Tstar)
|
|
# Using the leading constant from McLinden, 2000 since the leading term from Huber 2003 gives crazy values
|
|
eta_star = 26.692e-3 * sqrt(molemass * T) / (pow(sigma, 2) * OMEGA_2_2) / 1e6
|
|
return eta_star
|
|
|
|
|
|
def get_psi(fluid, ref_fluid, eta, T, rhomolar, e_k, sigma_nm):
|
|
|
|
THIS = CoolProp.AbstractState('HEOS', fluid)
|
|
REF = CoolProp.AbstractState('HEOS', ref_fluid)
|
|
|
|
THIS.update(CoolProp.DmolarT_INPUTS, rhomolar, T)
|
|
|
|
def residual_for_psi(psi, REF):
|
|
|
|
# Calculate the conformal state
|
|
conformal_state = THIS.conformal_state(ref_fluid, -1, -1)
|
|
# Calculate ESRR (which are based on the CONFORMAL state values)
|
|
f = THIS.T() / conformal_state['T'];
|
|
h = conformal_state['rhomolar'] / THIS.rhomolar(); # Must be the ratio of MOLAR densities!!
|
|
|
|
# The F factor
|
|
F_eta = sqrt(f) * pow(h, -2.0 / 3.0) * sqrt(THIS.molar_mass() / REF.molar_mass());
|
|
|
|
# Dilute viscosity of fluid of interest
|
|
eta_dilute = viscosity_dilute(fluid, T, e_k, sigma_nm)
|
|
|
|
# Required background contribution from reference fluid
|
|
viscosity_background_required = (eta - eta_dilute) / F_eta
|
|
|
|
REF.update(CoolProp.DmolarT_INPUTS, conformal_state['rhomolar'] * psi, conformal_state['T'])
|
|
visc_ref = REF.viscosity_contributions()
|
|
residual = visc_ref['initial_density'] + visc_ref['residual']
|
|
return residual - viscosity_background_required
|
|
|
|
psi = scipy.optimize.newton(residual_for_psi, 1.0, args=(REF,))
|
|
return psi
|
|
|
|
|
|
def arrayize(*args):
|
|
return [np.array(a) for a in args]
|
|
|
|
|
|
def HFO():
|
|
# Data from Zhao et al. dx.doi.org/10.1021/je5001457 | J. Chem. Eng. Data 2014, 59, 1366-1371
|
|
data_R1234yf = """293.15 1109.9 32.8 12.04 0.1442 6.82
|
|
303.09 1073.5 43.6 12.53 0.1319 5.71
|
|
313.20 1033.6 57.8 13.16 0.1223 4.60
|
|
323.19 990.2 76.0 13.88 0.1126 3.55
|
|
333.14 941.4 99.7 14.82 0.1016 2.55
|
|
343.11 883.5 132.2 16.12 0.0899 1.64
|
|
353.08 809.6 179.9 18.17 0.0820 0.81
|
|
358.05 761.5 214.8 19.78 0.0770 0.46
|
|
363.05 695.7 267.7 22.44 0.0700 0.15
|
|
365.05 657.4 301.0 24.26 0.0624 0.05"""
|
|
|
|
# Data from Zhao et al. dx.doi.org/10.1021/je5001457 | J. Chem. Eng. Data 2014, 59, 1366-1371
|
|
data_R1234zeE = """295.23 1172.5 24.1 12.11 0.1776 8.88
|
|
303.19 1146.1 30.6 12.46 0.1607 7.91
|
|
313.21 1111.1 40.8 12.93 0.1429 6.66
|
|
323.19 1073.6 53.6 13.46 0.1319 5.48
|
|
333.00 1033.3 69.8 14.06 0.1193 4.36
|
|
343.05 986.7 91.3 14.82 0.1132 3.30
|
|
353.00 924.0 119.7 15.80 0.1051 2.26
|
|
363.12 866.8 160.4 17.28 0.0924 1.35
|
|
373.14 776.9 225.2 19.89 0.0817 0.54"""
|
|
|
|
for fluid, data, e_k, sigma_nm in zip(['R1234yf', 'R1234ze(E)'], [data_R1234yf, data_R1234zeE], [281.14, 292.11], [0.5328, 0.5017]):
|
|
xx, yy, RHO, ETA, ETACP, ETARP = [], [], [], [], [], []
|
|
for line in data.split('\n'):
|
|
T, rhoL, rhoV, etaV, nuL, sigma = line.strip().split(' ')
|
|
rhoL = float(rhoL)
|
|
T = float(T)
|
|
nuL = float(nuL)
|
|
rhomolar = rhoL / CoolProp.CoolProp.PropsSI(fluid, 'molemass')
|
|
eta = nuL / 1000**2 * rhoL
|
|
psi = get_psi(fluid, 'Propane', eta, T, rhomolar, e_k, sigma_nm)
|
|
xx.append(T)
|
|
yy.append(psi)
|
|
RHO.append(rhomolar)
|
|
ETA.append(eta)
|
|
ETACP.append(CoolProp.CoolProp.PropsSI('V', 'T', T, 'Q', 0, fluid))
|
|
ETARP.append(CoolProp.CoolProp.PropsSI('V', 'T', T, 'Q', 0, 'REFPROP::' + CoolProp.CoolProp.get_fluid_param_string(fluid, 'REFPROP_name')))
|
|
|
|
RHO, xx, ETA, ETACP, ETARP = arrayize(RHO, xx, ETA, ETACP, ETARP)
|
|
rhor = RHO / CoolProp.CoolProp.PropsSI(fluid, 'rhomolar_critical')
|
|
|
|
plt.plot(rhor, yy, 'o-', label='from experimental data')
|
|
p = np.polyfit(rhor, yy, 2)
|
|
print(p[::-1])
|
|
plt.plot(rhor, np.polyval(p, rhor), 'o-', label='from correlation')
|
|
plt.xlabel(r'$\rho_r$')
|
|
plt.ylabel('$\psi$')
|
|
plt.legend(loc='best')
|
|
plt.show()
|
|
|
|
plt.title(fluid)
|
|
plt.plot(xx, (ETACP / ETA - 1) * 100, '^', label='CoolProp')
|
|
plt.plot(xx, (ETARP / ETA - 1) * 100, 'o', label='REFPROP')
|
|
plt.xlabel('Temperature (K)')
|
|
plt.ylabel('$100\\times(\eta_{calc}/\eta_{exp}-1)$ (%)')
|
|
plt.legend(loc='best')
|
|
plt.savefig(fluid + '_deviation.pdf')
|
|
plt.show()
|
|
|
|
|
|
def pentanes():
|
|
# from doi 10.1021/je0202174 | J. Chem. Eng. Data 2003, 48, 1418-1421
|
|
# T (K), rhoL (kg/m^3), rhoV (kg/m^3), eta (mPa-s)
|
|
data_cyclopentane = """253.15 258.15 263.15 268.15 273.15 278.15 283.15 288.15 293.15 298.15 303.15 308.15 313.15 318.15 323.15 328.15 333.15 338.15 343.15 348.15 353.15
|
|
784.64 779.53 774.59 769.77 765.12 760.20 755.32 750.27 745.02 738.63 731.97 725.15 718.32 711.59 705.11 699.08 693.40 688.44 684.25 680.96 678.71
|
|
0.0881 0.1127 0.1443 0.1848 0.2368 0.3036 0.3894 0.4999 0.6421 0.8255 1.062 1.368 1.764 2.279 2.950 3.827 4.980 6.509 8.554 11.33 15.20
|
|
0.7268 0.6786 0.6347 0.5930 0.5567 0.5224 0.4922 0.4646 0.4382 0.4148 0.3923 0.3714 0.3521 0.3350 0.3190 0.3048 0.2912 0.2793 0.2690 0.2590 0.2502"""
|
|
|
|
# from doi 10.1021/je0202174 | J. Chem. Eng. Data 2003, 48, 1418-1421
|
|
# T (K), rhoL (kg/m^3), rhoV (kg/m^3), eta (mPa-s)
|
|
data_isopentane = """253.15 258.15 263.15 268.15 273.15 278.15 283.15 288.15 293.15 298.15 303.15 308.15 313.15 318.15 323.15 328.15 333.15 338.15 343.15 348.15 353.15
|
|
658.32 653.55 648.73 643.87 639.01 634.15 629.35 624.63 620.05 615.69 610.87 605.63 600.05 594.23 588.24 582.18 576.13 570.18 564.41 558.92 553.79
|
|
0.4655 0.5889 0.7372 0.9137 1.122 1.366 1.650 1.979 2.356 2.788 3.278 3.833 4.459 5.162 5.949 6.827 7.803 8.886 10.09 11.41 12.87
|
|
0.3893 0.3661 0.3439 0.3201 0.3023 0.2859 0.2703 0.2547 0.2399 0.2289 0.2144 0.2023 0.1910 0.1813 0.1724 0.1611 0.1543 0.1480 0.1411 0.1332 0.1287"""
|
|
|
|
fluid = ''
|
|
|
|
def undelimit(args, delim=''):
|
|
return [np.array([float(_) for _ in a.strip().split(delim)]) for a in args]
|
|
|
|
from CoolProp.CoolProp import PropsSI
|
|
for fluid, e_k, sigma_nm in zip(['CycloPentane', 'Isopentane'], [406.33, 341.06], [0.518, 0.56232]):
|
|
xx, yy, RHO, ETA, ETACP, ETARP = [], [], [], [], [], []
|
|
for _T, _rhoLmass, _rhoVmass, _eta_mPas in zip(*undelimit(data_cyclopentane.split('\n'), delim=' ')):
|
|
MM = PropsSI('molemass', fluid)
|
|
rhomolar = _rhoLmass / MM
|
|
eta = _eta_mPas / 1000
|
|
psi = get_psi(fluid, 'Propane', eta, _T, rhomolar, e_k, sigma_nm)
|
|
xx.append(_T)
|
|
yy.append(psi)
|
|
RHO.append(rhomolar)
|
|
try:
|
|
ETACP.append(CoolProp.CoolProp.PropsSI('V', 'T', _T, 'Q', 0, fluid))
|
|
except:
|
|
ETACP.append(np.nan)
|
|
ETA.append(eta)
|
|
ETARP.append(CoolProp.CoolProp.PropsSI('V', 'T', _T, 'Q', 0, 'REFPROP::' + CoolProp.CoolProp.get_fluid_param_string(fluid, 'REFPROP_name')))
|
|
xx, yy, ETACP, ETARP, RHO, = arrayize(xx, yy, ETACP, ETARP, RHO)
|
|
rhored = CoolProp.CoolProp.PropsSI(fluid, 'rhomolar_critical')
|
|
print('rhored', rhored)
|
|
rhor = np.array(RHO) / rhored
|
|
|
|
plt.title(fluid)
|
|
plt.plot(rhor, yy, 'o-', label='from experimental data')
|
|
p = np.polyfit(rhor, yy, 2)
|
|
print(p[::-1])
|
|
plt.plot(rhor, np.polyval(p, rhor), 'o-', label='from correlation')
|
|
plt.xlabel(r'$\rho_r$')
|
|
plt.ylabel('$\psi$')
|
|
plt.legend(loc='best')
|
|
plt.show()
|
|
|
|
plt.title(fluid)
|
|
plt.plot(xx, (ETACP / ETA - 1) * 100, '^', label='CoolProp')
|
|
plt.plot(xx, (ETARP / ETA - 1) * 100, 'o', label='REFPROP')
|
|
plt.xlabel('Temperature (K)')
|
|
plt.ylabel('$100\\times(\eta_{calc}/\eta_{exp}-1)$ (%)')
|
|
plt.legend(loc='best')
|
|
plt.savefig(fluid + '_deviation.pdf')
|
|
plt.show()
|
|
|
|
|
|
HFO()
|
|
pentanes()
|