diff --git a/CoolPropBibTeXLibrary.bib b/CoolPropBibTeXLibrary.bib index 8b6e02b9..c22c0e5e 100644 --- a/CoolPropBibTeXLibrary.bib +++ b/CoolPropBibTeXLibrary.bib @@ -156,6 +156,21 @@ timestamp = {2013.04.28} } +@ARTICLE{Bell-IECR-2014, + author = {Bell, Ian H. and Wronski, Jorrit and Quoilin, Sylvain and Lemort, + Vincent}, + title = {Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and + the Open-Source Thermophysical Property Library CoolProp}, + journal = {Industrial \& Engineering Chemistry Research}, + year = {2014}, + volume = {53}, + pages = {2498-2508}, + number = {6}, + doi = {10.1021/ie4033999}, + eprint = {http://pubs.acs.org/doi/pdf/10.1021/ie4033999}, + url = {http://pubs.acs.org/doi/abs/10.1021/ie4033999} +} + @ARTICLE{Buecker-JPCRD-2006, author = {D. Buecker and W. Wagner}, title = {{A Reference Equation of State for the Thermodynamic Properties of diff --git a/dev/fluids/R124.json b/dev/fluids/R124.json index 780b4aa1..adae1115 100644 --- a/dev/fluids/R124.json +++ b/dev/fluids/R124.json @@ -279,5 +279,28 @@ "rhoVtriple_units": "mol/m^3" } ], - "NAME": "R124" + "NAME": "R124", + "TRANSPORT": { + "viscosity": { + "BibTeX": "", + "ECS_psi": { + "a": [ + 1.04253, + 0.00138528 + ], + "rhomolar_reducing": 4103.279546177282, + "rhomolar_reducing_units": "mol/m^3", + "t": [ + 0, + 1 + ] + }, + "ECS_reference_fluid": "Propane", + "epsilon_over_k": 275.8, + "epsilon_over_k_units": "K", + "sigma": 5.501e-10, + "sigma_units": "m", + "type": "ECS" + } + } } \ No newline at end of file diff --git a/include/AbstractState.h b/include/AbstractState.h index 981c5296..a090d783 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -221,6 +221,13 @@ protected: /// Using this backend, get the triple point temperature in K virtual long double calc_Ttriple(void){throw NotImplementedError("calc_Ttriple is not implemented for this backend");}; + /// Using this backend, get the critical point temperature in K + virtual long double calc_T_critical(void){throw NotImplementedError("calc_T_critical is not implemented for this backend");}; + /// Using this backend, get the critical point pressure in Pa + virtual long double calc_p_critical(void){throw NotImplementedError("calc_p_critical is not implemented for this backend");}; + /// Using this backend, get the critical point molar density in mol/m^3 + virtual long double calc_rhomolar_critical(void){throw NotImplementedError("calc_rhomolar_critical is not implemented for this backend");}; + public: virtual long double calc_melt_p_T(long double T){throw NotImplementedError("calc_melt_p_T is not implemented for this backend");}; @@ -271,6 +278,10 @@ public: double pmax(void); double Ttriple(void); + double T_critical(void); + double p_critical(void); + double rhomolar_critical(void); + std::string name(){return calc_name();}; // ---------------------------------------- @@ -402,7 +413,7 @@ public: /** -This class is a wrapper around an AbstractState. It handles construction and desctruction of the AbstractState +This class is a wrapper around an AbstractState. It handles construction and destruction of the AbstractState instance which decreases the likelihood of memory leaks. Only a few functions are exposed out of the AbstractState - this can be expanded, but functions must be manually exported out of the wrapper, which is not so nice. */ diff --git a/include/CoolPropFluid.h b/include/CoolPropFluid.h index 293457d2..44b0f0a0 100644 --- a/include/CoolPropFluid.h +++ b/include/CoolPropFluid.h @@ -189,6 +189,11 @@ struct ViscosityHigherOrderVariables ViscosityHigherOrderVariables(){type = VISCOSITY_HIGHER_ORDER_NOT_SET;}; }; +struct ViscosityECSVariables{ + long double rhomolar_reducing; + std::vector psi_a, psi_t; +}; + class TransportPropertyData { public: @@ -206,6 +211,7 @@ public: ViscosityDiluteVariables viscosity_dilute; ViscosityInitialDensityVariables viscosity_initial; ViscosityHigherOrderVariables viscosity_higher_order; + ViscosityECSVariables viscosity_ecs; ConductivityDiluteVariables conductivity_dilute; ConductivityResidualVariables conductivity_residual; ConductivityCriticalVariables conductivity_critical; diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 637a65b2..6ae8e244 100644 --- a/src/AbstractState.cpp +++ b/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(); diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.cpp b/src/Backends/Helmholtz/Fluids/FluidLibrary.cpp index c310aa85..9ddda4d3 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.cpp +++ b/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){ diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index a64a5541..fb006571 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/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 /** diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 5542e19e..30f1c688 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/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 { diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index b321000c..f734692d 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/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); diff --git a/src/Backends/Helmholtz/TransportRoutines.cpp b/src/Backends/Helmholtz/TransportRoutines.cpp index 55bc2c65..9ed3cbfa 100644 --- a/src/Backends/Helmholtz/TransportRoutines.cpp +++ b/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(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 &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(-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 */ \ No newline at end of file diff --git a/src/Backends/Helmholtz/TransportRoutines.h b/src/Backends/Helmholtz/TransportRoutines.h index fdecdfa2..3fe65fb9 100644 --- a/src/Backends/Helmholtz/TransportRoutines.h +++ b/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 */