diff --git a/dev/fluids/R124.json b/dev/fluids/R124.json index b67d1fc5..fd2fdb71 100644 --- a/dev/fluids/R124.json +++ b/dev/fluids/R124.json @@ -291,7 +291,9 @@ "TRANSPORT": { "viscosity": { "BibTeX": "", - "ECS_psi": { + "epsilon_over_k": 275.8, + "epsilon_over_k_units": "K", + "psi": { "a": [ 1.04253, 0.00138528 @@ -303,11 +305,9 @@ 1 ] }, - "ECS_reference_fluid": "Propane", - "epsilon_over_k": 275.8, - "epsilon_over_k_units": "K", - "sigma": 5.501e-10, - "sigma_units": "m", + "reference_fluid": "Propane", + "sigma_eta": 5.501e-10, + "sigma_eta_units": "m", "type": "ECS" } } diff --git a/include/CoolPropFluid.h b/include/CoolPropFluid.h index 44b0f0a0..9d8cba95 100644 --- a/include/CoolPropFluid.h +++ b/include/CoolPropFluid.h @@ -190,7 +190,8 @@ struct ViscosityHigherOrderVariables }; struct ViscosityECSVariables{ - long double rhomolar_reducing; + std::string reference_fluid; + long double psi_rhomolar_reducing; std::vector psi_a, psi_t; }; @@ -217,10 +218,13 @@ public: ConductivityCriticalVariables conductivity_critical; std::string BibTeX_viscosity, BibTeX_conductivity; + bool using_ECS; ///< A flag for whether to use extended corresponding states. False for no long double sigma_eta, epsilon_over_k; int hardcoded_viscosity, hardcoded_conductivity; TransportPropertyData(){hardcoded_viscosity = VISCOSITY_NOT_HARDCODED; - hardcoded_conductivity = CONDUCTIVITY_NOT_HARDCODED;}; + hardcoded_conductivity = CONDUCTIVITY_NOT_HARDCODED; + using_ECS = false; + }; }; /** @@ -513,6 +517,7 @@ class CoolPropFluid { EnvironmentalFactorsStruct environment; Ancillaries ancillaries; TransportPropertyData transport; + SimpleState crit; double gas_constant(){ return pEOS->R_u; }; double molar_mass(){ return pEOS->molar_mass; }; diff --git a/include/rapidjson/rapidjson_include.h b/include/rapidjson/rapidjson_include.h index 9d1cd3a6..ecf313a0 100644 --- a/include/rapidjson/rapidjson_include.h +++ b/include/rapidjson/rapidjson_include.h @@ -21,6 +21,14 @@ typedef unsigned int UINT32; namespace cpjson { + inline std::string json2string(rapidjson::Value &v) + { + rapidjson::StringBuffer buffer; + rapidjson::PrettyWriter writer(buffer); + + v.Accept(writer); + return buffer.GetString(); + } /// A convenience function to get a double from a JSON value, including error checking inline int get_integer(rapidjson::Value &v, std::string m) { diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index fb006571..a766ff45 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -437,11 +437,39 @@ protected: } }; + void parse_ECS_viscosity(rapidjson::Value &viscosity, CoolPropFluid & fluid) + { + fluid.transport.viscosity_ecs.reference_fluid = cpjson::get_string(viscosity,"reference_fluid"); + + // Parameters for correction polynomial + fluid.transport.viscosity_ecs.psi_a = cpjson::get_long_double_array(viscosity["psi"]["a"]); + fluid.transport.viscosity_ecs.psi_t = cpjson::get_long_double_array(viscosity["psi"]["t"]); + fluid.transport.viscosity_ecs.psi_rhomolar_reducing = cpjson::get_double(viscosity["psi"],"rhomolar_reducing"); + + fluid.transport.using_ECS = true; + } + /// Parse the transport properties void parse_viscosity(rapidjson::Value &viscosity, CoolPropFluid & fluid) { // Load the BibTeX key fluid.transport.BibTeX_viscosity = cpjson::get_string(viscosity,"BibTeX"); + + // Set the Lennard-Jones 12-6 potential variables, or approximate them from method of Chung + if (!viscosity.HasMember("sigma_eta")|| !viscosity.HasMember("epsilon_over_k")){ + default_transport(fluid); + } + else{ + fluid.transport.sigma_eta = cpjson::get_double(viscosity, "sigma_eta"); + fluid.transport.epsilon_over_k = cpjson::get_double(viscosity, "epsilon_over_k"); + } + + // If it is using ECS, set ECS parameters and quit + if (viscosity.HasMember("type") && !cpjson::get_string(viscosity, "type").compare("ECS")){ + parse_ECS_viscosity(viscosity, fluid); + return; + } + if (viscosity.HasMember("hardcoded")){ std::string target = cpjson::get_string(viscosity,"hardcoded"); if (!target.compare("Water")){ @@ -458,14 +486,7 @@ protected: } } - // Set the Lennard-Jones 12-6 potential variables, or approximate them from method of Chung - if (!viscosity.HasMember("sigma_eta")|| !viscosity.HasMember("epsilon_over_k")){ - default_transport(fluid); - } - else{ - fluid.transport.sigma_eta = cpjson::get_double(viscosity, "sigma_eta"); - fluid.transport.epsilon_over_k = cpjson::get_double(viscosity, "epsilon_over_k"); - } + // Load dilute viscosity term if (viscosity.HasMember("dilute")){ @@ -650,23 +671,6 @@ protected: /// Parse the transport properties void parse_transport(rapidjson::Value &transport, CoolPropFluid & fluid) { - //if (!fluid.name.compare("n-Hexane")){ - // rapidjson::Document dd; - // dd.SetObject(); - - // dd.AddMember("core",transport,dd.GetAllocator()); - // rapidjson::StringBuffer buffer; - // rapidjson::PrettyWriter writer(buffer); - - // dd.Accept(writer); - // std::string json0 = buffer.GetString(); - // //std::cout << json0 << std::endl; - - // FILE *fp; - // fp = fopen("nHexane_transport.json","w"); - // fprintf(fp,"%s",json0.c_str()); - // fclose(fp); - //} // Parse viscosity if (transport.HasMember("viscosity")){ @@ -691,8 +695,11 @@ protected: } /// Parse the critical state for the given EOS - void parse_crit_state(rapidjson::Value &alphar) + void parse_crit_state(rapidjson::Value &crit, CoolPropFluid & fluid) { + fluid.crit.T = cpjson::get_double(crit,"T"); + fluid.crit.p = cpjson::get_double(crit,"p"); + fluid.crit.rhomolar = cpjson::get_double(crit,"rhomolar"); }; /// Parse the critical state for the given EOS @@ -751,7 +758,8 @@ public: fluid.name = fluid_json["NAME"].GetString(); name_vector.push_back(fluid.name); // CAS number fluid.CAS = fluid_json["CAS"].GetString(); - + // Critical state + parse_crit_state(fluid_json["CRITICAL"], fluid); 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); @@ -761,7 +769,7 @@ public: fluid.aliases = cpjson::get_string_array(fluid_json["ALIASES"]); // EOS - parse_EOS_listing(fluid_json["EOS"],fluid); + parse_EOS_listing(fluid_json["EOS"], fluid); // Validate the fluid validate(fluid); @@ -789,9 +797,7 @@ 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); diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 30f1c688..4dcf839d 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -195,9 +195,18 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void) { if (is_pure_or_pseudopure) { - if (components[0]->transport.hardcoded_viscosity != CoolProp::TransportPropertyData::VISCOSITY_NOT_HARDCODED) + // Get a reference for code cleanness + CoolPropFluid &component = *(components[0]); + + // Check if using ECS + if (component.transport.using_ECS) { - switch(components[0]->transport.hardcoded_viscosity) + return TransportRoutines::viscosity_ECS(*this, *this); + } + + if (component.transport.hardcoded_viscosity != CoolProp::TransportPropertyData::VISCOSITY_NOT_HARDCODED) + { + switch(component.transport.hardcoded_viscosity) { case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_WATER: return TransportRoutines::viscosity_water_hardcoded(*this); @@ -206,7 +215,7 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void) case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_R23: return TransportRoutines::viscosity_R23_hardcoded(*this); default: - throw ValueError(format("hardcoded viscosity type [%d] is invalid for fluid %s", components[0]->transport.hardcoded_viscosity, name().c_str())); + throw ValueError(format("hardcoded viscosity type [%d] is invalid for fluid %s", component.transport.hardcoded_viscosity, name().c_str())); } } // Dilute part @@ -215,8 +224,9 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void) // 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 + // Critical part (no fluids have critical enhancement for viscosity currently) long double eta_critical = 0; + return eta_dilute + eta_back + eta_critical; } else @@ -311,6 +321,33 @@ std::string HelmholtzEOSMixtureBackend::calc_name(void) return components[0]->name; } } +long double HelmholtzEOSMixtureBackend::calc_T_critical(void) +{ + if (components.size() != 1){ + throw ValueError(format("For now, calc_T_critical is only valid for pure and pseudo-pure fluids, %d components", components.size())); + } + else{ + return components[0]->crit.T; + } +} +long double HelmholtzEOSMixtureBackend::calc_p_critical(void) +{ + if (components.size() != 1){ + throw ValueError(format("For now, calc_p_critical is only valid for pure and pseudo-pure fluids, %d components", components.size())); + } + else{ + return components[0]->crit.p; + } +} +long double HelmholtzEOSMixtureBackend::calc_rhomolar_critical(void) +{ + if (components.size() != 1){ + throw ValueError(format("For now, calc_rhomolar_critical is only valid for pure and pseudo-pure fluids, %d components", components.size())); + } + else{ + return components[0]->crit.rhomolar; + } +} long double HelmholtzEOSMixtureBackend::calc_Tmax(void) { double summer = 0; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index f734692d..d21339eb 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -128,6 +128,10 @@ public: long double calc_pmax(void); long double calc_Ttriple(void); + long double calc_T_critical(void); + long double calc_p_critical(void); + long double calc_rhomolar_critical(void); + std::string calc_name(void); long double calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector & mole_fractions, const long double &tau, const long double &delta); diff --git a/src/Backends/Helmholtz/TransportRoutines.cpp b/src/Backends/Helmholtz/TransportRoutines.cpp index 9ed3cbfa..3b76b6b5 100644 --- a/src/Backends/Helmholtz/TransportRoutines.cpp +++ b/src/Backends/Helmholtz/TransportRoutines.cpp @@ -831,10 +831,10 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H std::vector &a = ECS.psi_a, &t = ECS.psi_t; - // The correction polynomial + // The correction polynomial psi_eta double psi = 0; - for (std::size_t i=0; i<=a.size(); i++) - psi += a[i]*pow(HEOS.rhomolar()/ECS.rhomolar_reducing, t[i]); + for (std::size_t i=0; i < a.size(); i++) + psi += a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, t[i]); // The dilute gas portion for the fluid of interest [Pa-s] long double eta_dilute = viscosity_dilute_kinetic_theory(HEOS);