Baby steps towards transport properties.

src/Backends/Helmholtz/FlashRoutines.h

@@ -41,17 +41,17 @@ public:
     
     /// A generic flash routine for the pairs (T,D), (T,H), (T,S), and (T,U).  Similar analysis is needed
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
-    /// @param other The index for the other input, see CoolProp::parameters; allowed values are iDmolar, iHmolar, iSmolar, iUmolar
+    /// @param other The index for the other input from CoolProp::parameters; allowed values are iDmolar, iHmolar, iSmolar, iUmolar
     static void DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, int other);
     
     /// A generic flash routine for the pairs (P,H), (P,S), and (P,U).  Similar analysis is needed
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
-    /// @param other The index for the other input, see CoolProp::parameters; allowed values are iHmolar, iSmolar, iUmolar
+    /// @param other The index for the other input from CoolProp::parameters; allowed values are iHmolar, iSmolar, iUmolar
     static void HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other);
     
     /// A generic flash routine for the pairs (D,P), (D,H), (D,S), and (D,U).  Similar analysis is needed
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
-    /// @param other The index for the other input, see CoolProp::parameters; allowed values are iP, iHmolar, iSmolar, iUmolar
+    /// @param other The index for the other input from CoolProp::parameters; allowed values are iP, iHmolar, iSmolar, iUmolar
     static void PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other);
 };
 

src/Backends/Helmholtz/Fluids/CoolPropFluid.h

@@ -58,6 +58,13 @@ public:
     double residual(double T, double rhomolar);
     double critical(double T, double rhomolar);
 };
+class TransportPropertyData
+{
+public:
+    ViscosityCorrelation viscosity;
+    ThermalConductivityCorrelation conductivity;
+    long double sigma_eta, epsilon_over_k;
+};
 
 /**
 The surface tension correlation class uses correlations for the surface tension that are all
@@ -184,6 +191,7 @@ public:
     }
 };
 
+
 class MeltingLine
 {
 };
@@ -343,6 +351,7 @@ class CoolPropFluid {
         BibTeXKeysStruct BibTeXKeys;
         EnvironmentalFactorsStruct environment;
         Ancillaries ancillaries;
+        TransportPropertyData transport;
 
         double gas_constant(){ return pEOS->R_u; };
         double molar_mass(){ return pEOS->molar_mass; };

src/Backends/Helmholtz/Fluids/FluidLibrary.h

@@ -257,9 +257,22 @@ protected:
         fluid.pEOS = &(fluid.EOSVector[0]);
     };
 
-    /// Parse the reducing state for the given EOS
-    void parse_reducing_state(rapidjson::Value &alphar)
+    /// Parse the transport properties
+    void parse_transport(rapidjson::Value &transport, CoolPropFluid & fluid)
     {
+        if (!transport.HasMember("sigma_eta")|| !transport.HasMember("epsilon_over_k")){
+            // Use the method of Chung to approximate the values for epsilon_over_k and sigma_eta
+            // Chung, T.-H.; Ajlan, M.; Lee, L. L.; Starling, K. E. Generalized Multiparameter Correlation for Nonpolar and Polar Fluid Transport Properties. Ind. Eng. Chem. Res. 1988, 27, 671-679.
+            // rhoc needs to be in mol/L to yield a sigma in nm, 
+            long double rho_crit_molar = fluid.pEOS->reduce.rhomolar/1000.0;// [mol/m3 to mol/L]
+            long double Tc = fluid.pEOS->reduce.T;
+            fluid.transport.sigma_eta = 0.809/pow(rho_crit_molar, static_cast<long double>(1.0/3.0))/1e9; // 1e9 is to convert from nm to m
+            fluid.transport.epsilon_over_k = Tc/1.3593; // [K]
+        }
+        else{
+            fluid.transport.sigma_eta = cpjson::get_double(transport, "sigma_eta");
+            fluid.transport.epsilon_over_k = cpjson::get_double(transport, "epsilon_over_k");
+        }
     };
 
     /// Parse the critical state for the given EOS
@@ -352,6 +365,14 @@ public:
         else{
             parse_environmental(fluid_json["ENVIRONMENTAL"], fluid);
         }
+
+        // Parse the environmental parameters
+        if (!(fluid_json.HasMember("TRANSPORT"))){
+            std::cout << format("Transport property data are missing for fluid [%s]\n", fluid.name.c_str()) ;
+        }
+        else{
+            parse_transport(fluid_json["TRANSPORT"], fluid);
+        }
         
         // If the fluid is ok...
 

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -36,6 +36,7 @@ public:
     ExcessTerm Excess;
 
     friend class FlashRoutines; // Allows the routines in the FlashRoutines class to have access to all the protected members and methods of this class
+    friend class TransportRoutines; // Allows the routines in the TransportRoutines class to have access to all the protected members and methods of this class
 
     // Helmholtz EOS backend uses mole fractions
     bool using_mole_fractions(){return true;}

src/Backends/Helmholtz/TransportRoutines.cpp

@@ -0,0 +1,27 @@
+
+#include "TransportRoutines.h"
+
+namespace CoolProp{
+
+long double TransportRoutines::general_dilute_gas_viscosity(HelmholtzEOSMixtureBackend &HEOS)
+{
+    if (HEOS.is_pure_or_pseudopure)
+    {
+        long double Tstar = HEOS.T()/HEOS.components[0]->transport.epsilon_over_k;
+        long double sigma_nm = HEOS.components[0]->transport.sigma_eta*1e9; // 1e9 to convert from m to nm
+        long double molar_mass_kgkmol = HEOS.molar_mass()*1000; // 1000 to convert from kg/mol to kg/kmol
+    
+        // The nondimensional empirical collision integral from Neufeld 
+        // Neufeld, P. D.; Janzen, A. R.; Aziz, R. A. Empirical Equations to Calculate 16 of the Transport Collision Integrals (l,s)* 
+        // for the Lennard-Jones (12-6) Potential. J. Chem. Phys. 1972, 57, 1100-1102
+        long double OMEGA22 = 1.16145*pow(Tstar, static_cast<long double>(-0.14874))+0.52487*exp(-0.77320*Tstar)+2.16178*exp(-2.43787*Tstar);
+
+        // The dilute gas component - 
+        return 26.692e-9*sqrt(molar_mass_kgkmol*HEOS.T())/(pow(sigma_nm, 2)*OMEGA22); // Pa-s
+    }
+    else{
+        throw NotImplementedError("TransportRoutines::GeneralDiluteGasViscosity is only for pure and pseudo-pure");
+    }
+}
+
+}; /* namespace CoolProp */

src/Backends/Helmholtz/TransportRoutines.h

@@ -0,0 +1,39 @@
+#ifndef TRANSPORTROUTINES_H
+#define TRANSPORTROUTINES_H
+
+#include "HelmholtzEOSMixtureBackend.h"
+
+namespace CoolProp{
+
+class TransportRoutines
+{
+public:
+    /**
+    \brief The general dilute gas viscosity from used for ECS
+
+    \f[
+    \eta^0 = \displaystyle\frac{26.692\times 10^{-9}\sqrt{MT}}{\sigma^2\Omega^{(2,2)}(T^*)}
+    \f]
+    \f[
+    \Omega^{(2,2)}(T^*)=1.16145(T^*)^{-0.14874}+0.52487\exp(-0.77320T^*)+2.16178\exp(-2.43787T^*)
+    \f]
+    with \f$T^* = \frac{T}{\varepsilon/k}\f$ and \f$\sigma\f$ in nm, M is in kg/kmol. Yields viscosity in Pa-s.
+    */
+    static long double general_dilute_gas_viscosity(HelmholtzEOSMixtureBackend &HEOS);
+    
+    /**
+    \brief A dilute gas term that has a form like
+
+    \f[
+    \eta^0 = \displaystyle\frac{A\sqrt{MT}}{\sigma^2\mathfrak{S}(T^*)}
+    \f]
+    \f[
+    \mathfrak{S}(T^*)=\exp\left(\sum_ia_i[\ln T^*]^i\right)
+    \f]
+    with \f$T^* = \frac{T}{\varepsilon/k}\f$ and \f$\sigma\f$ in nm, M is in kg/kmol. Yields viscosity in Pa-s.
+    */
+    static long double dilute_gas_viscosity_collision_integral(HelmholtzEOSMixtureBackend &HEOS);
+};
+
+}; /* namespace CoolProp */
+#endif