First cut at adding code for ECS for viscosity. Starting with R124

src/AbstractState.cpp

@@ -242,6 +242,15 @@ double AbstractState::Ttriple(void){
 double AbstractState::pmax(void){
     return calc_pmax();
 }
+double AbstractState::T_critical(void){
+    return calc_T_critical();
+}
+double AbstractState::p_critical(void){
+    return calc_p_critical();
+}
+double AbstractState::rhomolar_critical(void){
+    return calc_rhomolar_critical();
+}
 
 double AbstractState::hmolar(void){
     if (!_hmolar) _hmolar = calc_hmolar();

src/Backends/Helmholtz/Fluids/FluidLibrary.cpp

@@ -11,7 +11,11 @@ void load()
 	rapidjson::Document dd;
     // This json formatted string comes from the all_fluids_JSON.h header which is a C++-escaped version of the JSON file
     dd.Parse<0>(all_fluids_JSON.c_str());
-	if (dd.HasParseError()){throw ValueError("Unable to load all_fluids.json");} else{library.add_many(dd);}
+	if (dd.HasParseError()){
+        throw ValueError("Unable to load all_fluids.json");
+    } else{
+        try{library.add_many(dd);}catch(std::exception &e){std::cout << e.what() << std::endl;}
+    }
 }
 
 JSONFluidLibrary & get_library(void){

src/Backends/Helmholtz/Fluids/FluidLibrary.h

@@ -752,6 +752,11 @@ public:
         // CAS number
         fluid.CAS = fluid_json["CAS"].GetString();
 
+
+        if (get_debug_level() > 5){
+            std::cout << format("Loading fluid %s with CAS %s; %d fluids loaded\n", fluid.name.c_str(), fluid.CAS.c_str(), index);
+        }
+
         // Aliases
         fluid.aliases = cpjson::get_string_array(fluid_json["ALIASES"]);
         
@@ -784,7 +789,9 @@ public:
         else{
             parse_environmental(fluid_json["ENVIRONMENTAL"], fluid);
         }
-
+        if (index == 80){
+            double rr =0;
+        }
         // Parse the environmental parameters
         if (!(fluid_json.HasMember("TRANSPORT"))){
             default_transport(fluid);
@@ -807,6 +814,7 @@ public:
             string_to_index_map[fluid.aliases[i]] = index;
         }
 
+        if (get_debug_level() > 5){ std::cout << format("Loaded.\n"); }
     };
     /// Get a CoolPropFluid instance stored in this library
     /**

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

@@ -156,6 +156,41 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity_dilute(void)
     }
     
 }
+long double HelmholtzEOSMixtureBackend::calc_viscosity_background()
+{
+    long double eta_dilute = calc_viscosity_dilute();
+    return calc_viscosity_background(eta_dilute);
+}
+long double HelmholtzEOSMixtureBackend::calc_viscosity_background(long double eta_dilute)
+{
+    // Residual part
+    long double B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this);
+    long double rho = rhomolar();
+    long double initial_part = eta_dilute*B_eta_initial*rhomolar();
+
+    // Higher order terms
+    long double delta_eta_h;
+    switch(components[0]->transport.viscosity_higher_order.type)
+    {
+    case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BATSCHINKI_HILDEBRAND:
+        delta_eta_h = TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand(*this); break;
+    case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_FRICTION_THEORY:
+        delta_eta_h = TransportRoutines::viscosity_higher_order_friction_theory(*this); break;
+    case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HYDROGEN:
+        delta_eta_h = TransportRoutines::viscosity_hydrogen_higher_order_hardcoded(*this); break;
+    case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HEXANE:
+        delta_eta_h = TransportRoutines::viscosity_hexane_higher_order_hardcoded(*this); break;
+    case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_ETHANE:
+        delta_eta_h = TransportRoutines::viscosity_ethane_higher_order_hardcoded(*this); break;
+    default:
+        throw ValueError(format("higher order viscosity type [%d] is invalid for fluid %s", components[0]->transport.viscosity_dilute.type, name().c_str()));
+    }
+
+    long double eta_residual = initial_part + delta_eta_h;
+
+    return eta_residual;
+}
+
 long double HelmholtzEOSMixtureBackend::calc_viscosity(void)
 {
     if (is_pure_or_pseudopure)
@@ -176,35 +211,13 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void)
         }
         // Dilute part
         long double eta_dilute = calc_viscosity_dilute();
-        
-        // Residual part
-        long double B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this);
-        long double rho = rhomolar();
-        long double initial_part = eta_dilute*B_eta_initial*rhomolar();
 
-        // Higher order terms
-        long double delta_eta_h;
-        switch(components[0]->transport.viscosity_higher_order.type)
-        {
-        case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BATSCHINKI_HILDEBRAND:
-            delta_eta_h = TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand(*this); break;
-        case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_FRICTION_THEORY:
-            delta_eta_h = TransportRoutines::viscosity_higher_order_friction_theory(*this); break;
-        case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HYDROGEN:
-            delta_eta_h = TransportRoutines::viscosity_hydrogen_higher_order_hardcoded(*this); break;
-        case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HEXANE:
-            delta_eta_h = TransportRoutines::viscosity_hexane_higher_order_hardcoded(*this); break;
-        case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_ETHANE:
-            delta_eta_h = TransportRoutines::viscosity_ethane_higher_order_hardcoded(*this); break;
-        default:
-            throw ValueError(format("higher order viscosity type [%d] is invalid for fluid %s", components[0]->transport.viscosity_dilute.type, name().c_str()));
-        }
-
-        long double eta_residual = initial_part + delta_eta_h;
+        // Background viscosity given by the sum of the initial density dependence and higher order terms
+        long double eta_back = calc_viscosity_background(eta_dilute);
 
         // Critical part
         long double eta_critical = 0;
-        return eta_dilute + eta_residual + eta_critical;
+        return eta_dilute + eta_back + eta_critical;
     }
     else
     {

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -35,8 +35,9 @@ public:
     ReducingFunctionContainer Reducing;
     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
+    friend class FlashRoutines; // Allows the static methods in the FlashRoutines class to have access to all the protected members and methods of this class
+    friend class TransportRoutines; // Allows the static methods in the TransportRoutines class to have access to all the protected members and methods of this class
+    //friend class MixtureDerivatives; // 
 
     // Helmholtz EOS backend uses mole fractions
     bool using_mole_fractions(){return true;}
@@ -119,6 +120,8 @@ public:
     long double calc_surface_tension(void);
     long double calc_viscosity(void);
     long double calc_viscosity_dilute(void);
+    long double calc_viscosity_background(void);
+    long double calc_viscosity_background(long double eta_dilute);
     long double calc_conductivity(void);
 
     long double calc_Tmax(void);

src/Backends/Helmholtz/TransportRoutines.cpp

@@ -584,6 +584,10 @@ long double TransportRoutines::conductivity_dilute_hardcoded_CO2(HelmholtzEOSMix
 	//Vesovic Eq. 12 [no units]
 	double r = sqrt(2.0/5.0*cint_k);
 
+    // According to REFPROP, 1+r^2 = cp-2.5R.  This is unclear to me but seems to suggest that cint/k is the difference 
+    // between the ideal gas specific heat and a monatomic specific heat of 5/2*R. Using the form of cint/k from Vesovic 
+    // does not yield exactly the correct values
+
 	Tstar = HEOS.T()/e_k;
 	//Vesovic Eq. 30 [no units]
 	summer = 0;
@@ -615,7 +619,7 @@ long double TransportRoutines::conductivity_dilute_eta0_and_poly(HelmholtzEOSMix
 
         double eta0_uPas = HEOS.calc_viscosity_dilute()*1e6; // [uPa-s]
         double summer = E.A[0]*eta0_uPas;
-	    for (int i=1; i < E.A.size(); ++i)
+	    for (std::size_t i=1; i < E.A.size(); ++i)
 		    summer += E.A[i]*pow(static_cast<long double>(HEOS.tau()), E.t[i]);
 	    return summer;
     }
@@ -632,7 +636,7 @@ long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureB
 				{-1.21051378,1.60812989,-0.621178141,0.0716373224,0,0},
 				{-2.7203370,4.57586331,-3.18369245,1.1168348,-0.19268305,0.012913842}};
 
-	double lambdabar_0,lambdabar_1,lambdabar_2,rhobar,Tbar,sum,R_Water;
+	double lambdabar_0,lambdabar_1,lambdabar_2,rhobar,Tbar,sum;
 	double Tstar=647.096,rhostar=322,pstar=22064000,lambdastar=1e-3,mustar=1e-6;
 	double tau,xi;
 	int i,j;
@@ -770,7 +774,7 @@ long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixture
     // Critical component
     lambda_c = 0.0;
     
-    if (3.5 < T & T < 12)
+    if (3.5 < T && T < 12)
     {
         double x0 = 0.392, E1 = 2.8461, E2 = 0.27156, beta = 0.3554, gamma = 1.1743, delta = 4.304, rhoc_crit = 69.158, 
             Tc = 5.18992, pc = 2.2746e5, R = 4.633e-10, m = 6.6455255e-27, k = 1.38066e-23, pi = M_PI;
@@ -813,4 +817,56 @@ long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixture
 	return lambda_0+lambda_e+lambda_c;
 }
 
+long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference)
+{
+    // Collect some parameters
+    long double M = HEOS.molar_mass(),
+                M0 = HEOS_Reference.molar_mass(),
+                Tc = HEOS.T_critical(),
+                Tc0 = HEOS_Reference.T_critical(),
+                rhocmolar = HEOS.rhomolar_critical(),
+                rhocmolar0 = HEOS_Reference.rhomolar_critical();
+
+    CoolProp::ViscosityECSVariables &ECS = HEOS.components[0]->transport.viscosity_ecs;
+
+    std::vector<long double> &a = ECS.psi_a, &t = ECS.psi_t;
+    
+    // The correction polynomial
+    double psi = 0;
+    for (std::size_t i=0; i<=a.size(); i++)
+        psi += a[i]*pow(HEOS.rhomolar()/ECS.rhomolar_reducing, t[i]);
+
+    // The dilute gas portion for the fluid of interest [Pa-s]
+    long double eta_dilute = viscosity_dilute_kinetic_theory(HEOS);
+
+    // Calculate the correction polynomial
+
+
+    /// \todo To be solved for...
+    // TODO: To be solved for...
+    long double theta = 1;
+    long double phi = 1;
+
+    // The equivalent substance reducing ratios
+    long double f = Tc/Tc0*theta;
+	long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
+
+    // To be solved for
+    long double T0 = HEOS.T()/f;
+    long double rhomolar0 = HEOS.rhomolar()*h;
+
+    // Update the reference fluid with the conformal state
+    HEOS_Reference.update(DmolarT_INPUTS, rhomolar0*psi, T0);
+
+    // The reference fluid's contribution to the viscosity [Pa-s]
+    long double eta_resid = HEOS_Reference.calc_viscosity_background();
+
+    // The F factor
+	long double F_eta = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M/M0);
+
+    // The total viscosity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
+	long double eta = eta_dilute + eta_resid*F_eta;
+
+    return eta;
+}
 }; /* namespace CoolProp */

src/Backends/Helmholtz/TransportRoutines.h

@@ -140,6 +140,8 @@ public:
     \f$\rho\f$ and \f$\rho_c\f$ are in mol\f$\cdot\f$m\f$^{-3}\f$, \f$\eta\f$ is the viscosity in Pa\f$\cdot\f$s, 
     and the remaining parameters are defined in the following tables.
 
+    It should be noted that some authors use slightly different values for the "universal" constants
+
     Coefficients for use in the simplified Olchowy-Sengers critical term
     Parameter             | Variable     | Value
     ---------             | --------     | ------
@@ -172,6 +174,25 @@ public:
 
     static long double conductivity_critical_hardcoded_ammonia(HelmholtzEOSMixtureBackend &HEOS);
 
+    /**
+    \brief Calculate the viscosity using the extended corresponding states method
+
+    This method is covered in depth in
+
+    Bell, I. H.; Wronski, J.; Quoilin, S. & Lemort, V. (2014), Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp, Industrial & Engineering Chemistry Research, 53, (6), 2498-2508
+
+    which is originally based on the methods presented in 
+
+    Huber, M. L., Laesecke, A. and Perkins, R. A., (2003), Model for the Viscosity and Thermal Conductivity of Refrigerants, Including a New Correlation for the Viscosity of R134a, Industrial & Engineering Chemistry Research, v. 42, pp. 3163-3178
+
+    and
+
+    McLinden, M. O.; Klein, S. A. & Perkins, R. A. (2000), An extended corresponding states model for the thermal conductivity of refrigerants and refrigerant mixtures, Int. J. Refrig., 23, 43-63
+
+
+    */
+    static long double viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference);
+
 
 }; /* class TransportRoutines */