Add viscosity fits for cyclo and isopentane; see #786

CoolProp · Nov 14, 2015 · 858714b · 858714b
1 parent 7151406
commit 858714b
Show file tree

Hide file tree

Showing 4 changed files with 183 additions and 64 deletions.
diff --git a/dev/fluids/Cyclopentane.json b/dev/fluids/Cyclopentane.json
@@ -448,5 +448,31 @@
       "smolar": 124.2473875963392, 
       "smolar_units": "J/mol/K"
     }
+  }, 
+  "TRANSPORT": {
+    "viscosity": {
+      "BibTeX": "BELL-PERSONAL-2015", 
+      "_note": "fit based on the values in Ma et al., JCED, 2003, doi 10.1021/je0202174 ", 
+      "epsilon_over_k": 406.33, 
+      "epsilon_over_k_units": "K", 
+      "psi": {
+        "a": [
+          1.61077063, 
+          -0.36230585, 
+          0.05748161
+        ], 
+        "rhomolar_reducing": 3820, 
+        "rhomolar_reducing_units": "mol/m^3", 
+        "t": [
+          0, 
+          1, 
+          2
+        ]
+      }, 
+      "reference_fluid": "n-Propane", 
+      "sigma_eta": 5.18e-10, 
+      "sigma_eta_units": "m", 
+      "type": "ECS"
+    }
   }
 }
diff --git a/dev/fluids/Isopentane.json b/dev/fluids/Isopentane.json
@@ -426,5 +426,31 @@
       "smolar": 176.18484251588657, 
       "smolar_units": "J/mol/K"
     }
+  }, 
+  "TRANSPORT": {
+    "viscosity": {
+      "BibTeX": "BELL-PERSONAL-2015", 
+      "_note": "fit based on the values in Ma et al., JCED, 2003, doi 10.1021/je0202174 ", 
+      "epsilon_over_k": 341.06, 
+      "epsilon_over_k_units": "K", 
+      "psi": {
+        "a": [
+          1.59193594, 
+          -0.36852042, 
+          0.04951376
+        ], 
+        "rhomolar_reducing": 3271.0, 
+        "rhomolar_reducing_units": "mol/m^3", 
+        "t": [
+          0, 
+          1, 
+          2
+        ]
+      }, 
+      "reference_fluid": "n-Propane", 
+      "sigma_eta": 5.6232e-10, 
+      "sigma_eta_units": "m", 
+      "type": "ECS"
+    }
   }
 }
diff --git a/dev/fluids/R1234ze(E).json b/dev/fluids/R1234ze(E).json
@@ -743,9 +743,9 @@
       "epsilon_over_k_units": "K", 
       "psi": {
         "a": [
-          0.77371022, 
-          0.16492734, 
-          -0.02851852
+          0.72309634, 
+          0.2235595, 
+          -0.04217
         ], 
         "rhomolar_reducing": 4290, 
         "rhomolar_reducing_units": "mol/m^3", 

diff --git a/dev/scripts/ECS_fitter.py b/dev/scripts/ECS_fitter.py
@@ -23,10 +23,8 @@ def get_psi(fluid, ref_fluid, eta, T, rhomolar, e_k, sigma_nm):
 
     THIS = CoolProp.AbstractState('HEOS', fluid)
     REF = CoolProp.AbstractState('HEOS', ref_fluid)
-    RP = CoolProp.AbstractState('REFPROP', fluid)
 
     THIS.update(CoolProp.DmolarT_INPUTS, rhomolar, T)
-    RP.update(CoolProp.DmolarT_INPUTS, rhomolar, T)
 
     def residual_for_psi(psi, REF):
 
@@ -53,63 +51,132 @@ def residual_for_psi(psi, REF):
     psi = scipy.optimize.newton(residual_for_psi, 1.0, args = (REF,))
     return psi
 
-# 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"""
-
 def arrayize(*args):
-    return [np.array(a) for a in args]
-
-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('R1234yf', '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-')
-    p = np.polyfit(rhor, yy, 2)
-    print p[::-1]
-    plt.plot(rhor, np.polyval(p, rhor), 'o-')
-    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()
+        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()