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