From f93aa762103ec17186456c9ed7ca2cd0f65496ca Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 28 May 2014 19:15:34 +0200 Subject: [PATCH] ECS conductivity work - nearly works, just need to sort out the conformal state solver. Goodbye to AbstractStateWrapper - can use std::tr1::shared_ptr, much nicer Signed-off-by: Ian Bell --- dev/README.md | 12 +- dev/fluids/R124.json | 29 ++++ include/AbstractState.h | 29 ---- include/CoolPropFluid.h | 15 +- src/AbstractState.cpp | 23 +-- src/Backends/Helmholtz/Fluids/FluidLibrary.h | 23 ++- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 75 ++++++--- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 2 + src/Backends/Helmholtz/TransportRoutines.cpp | 85 +++++++++- src/Backends/Helmholtz/TransportRoutines.h | 2 + src/HumidAirProp.cpp | 151 +++++++++--------- 11 files changed, 287 insertions(+), 159 deletions(-) diff --git a/dev/README.md b/dev/README.md index 4db7a04a..aedde6b0 100644 --- a/dev/README.md +++ b/dev/README.md @@ -30,4 +30,14 @@ JSON fluid definitions - the second make run will actually generate the program. Testing ------- -CMake generates a target for testing. You can build the test executable with `make testRunner`. \ No newline at end of file +CMake generates a target for testing. You can build the test executable with `make testRunner`. + +Geting Boost shared_ptr +----------------------- + +Download Boost sources +Expand zip file +Command prompt in root + + bootstrap + b2 tools/bcp diff --git a/dev/fluids/R124.json b/dev/fluids/R124.json index fd2fdb71..83dbb172 100644 --- a/dev/fluids/R124.json +++ b/dev/fluids/R124.json @@ -289,6 +289,35 @@ ], "NAME": "R124", "TRANSPORT": { + "conductivity": { + "BibTeX": "", + "f_int": { + "T_reducing": 1.0, + "T_reducing_units": "K", + "a": [ + 0.0011769, + 6.78397e-07 + ], + "t": [ + 0, + 1 + ] + }, + "psi": { + "a": [ + 1.0898, + -0.0154229 + ], + "rhomolar_reducing": 4103.279546177282, + "rhomolar_reducing_units": "mol/m^3", + "t": [ + 0, + 1 + ] + }, + "reference_fluid": "Propane", + "type": "ECS" + }, "viscosity": { "BibTeX": "", "epsilon_over_k": 275.8, diff --git a/include/AbstractState.h b/include/AbstractState.h index a090d783..7e15a1dd 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -411,34 +411,5 @@ public: */ }; - -/** -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. -*/ -class AbstractStateWrapper -{ -protected: - AbstractState *p; - // Copying is disabled for now until we determine the right semantics for ownership of AbstractState instance - AbstractStateWrapper(const AbstractStateWrapper& copy_from_me){}; -public: - AbstractStateWrapper(){this->p = NULL;}; - void set(const std::string &backend, const std::string &fluid_string){ - this->p = AbstractState::factory(backend, fluid_string); - }; - ~AbstractStateWrapper(){delete this->p;}; - void update(long input_pair, double Value1, double Value2){ - if (this->p == NULL) { throw ValueError("AbstractState in AbstractStateWrapper has not been instantiated yet");} - this->p->update(input_pair,Value1,Value2); - }; - double keyed_output(int key) { - if (this->p == NULL) { throw ValueError("AbstractState in AbstractStateWrapper has not been instantiated yet");} - return this->p->keyed_output(key); - } - bool empty(){return (this->p == NULL);} -}; - } /* namespace CoolProp */ #endif /* ABSTRACTSTATE_H_ */ diff --git a/include/CoolPropFluid.h b/include/CoolPropFluid.h index 9d8cba95..7daad769 100644 --- a/include/CoolPropFluid.h +++ b/include/CoolPropFluid.h @@ -44,6 +44,12 @@ struct EOSLimits double Tmin, Tmax, rhomax, pmax; }; +struct ConductivityECSVariables{ + std::string reference_fluid; + long double psi_rhomolar_reducing, f_int_T_reducing; + std::vector psi_a, psi_t, f_int_a, f_int_t; +}; + struct ConductivityDiluteEta0AndPolyData{ std::vector A, t; }; @@ -102,7 +108,7 @@ struct ConductivityCriticalSimplifiedOlchowySengersData{ // Suggested default values - can be over-written GAMMA = 0.0496; //[-] zeta0 = 1.94e-10; //[m] - qD = 1e9; //[m] + qD = 2e9; //[m] // Set to invalid number, can be provided in the JSON file T_ref = _HUGE; @@ -216,14 +222,17 @@ public: ConductivityDiluteVariables conductivity_dilute; ConductivityResidualVariables conductivity_residual; ConductivityCriticalVariables conductivity_critical; + ConductivityECSVariables conductivity_ecs; std::string BibTeX_viscosity, BibTeX_conductivity; - bool using_ECS; ///< A flag for whether to use extended corresponding states. False for no + bool viscosity_using_ECS; ///< A flag for whether to use extended corresponding states for viscosity. False for no + bool conductivity_using_ECS; ///< A flag for whether to use extended corresponding states for conductivity. 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; - using_ECS = false; + viscosity_using_ECS = false; + conductivity_using_ECS = false; }; }; diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 6ae8e244..d99f9e67 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -332,28 +332,7 @@ TEST_CASE("Check AbstractState","[AbstractState]") { SECTION("bad backend") { - CoolProp::AbstractStateWrapper Water; - CHECK_THROWS(Water.set("DEFINITELY_A_BAD_BACKEND", "Water")); - } -} - -TEST_CASE("Check AbstractStateWrapper","[AbstractStateWrapper]") -{ - SECTION("empty on init") - { - CoolProp::AbstractStateWrapper Water; - CHECK_NOTHROW(Water.empty()); - CHECK(Water.empty() == true); - CHECK_THROWS(Water.update(CoolProp::QT_INPUTS,1,300)); - } - SECTION("initialized") - { - CoolProp::AbstractStateWrapper Water; - CHECK_NOTHROW(Water.empty()); - CHECK(Water.empty() == true); - Water.set("HEOS", "Water"); - CHECK(Water.empty() == false); - CHECK_NOTHROW(Water.update(CoolProp::QT_INPUTS,1,300)); + CHECK_THROWS(std::tr1::shared_ptr Water(CoolProp::AbstractState::factory("DEFINITELY_A_BAD_BACKEND", "Water"))); } } diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index a766ff45..7b66fe96 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -437,6 +437,21 @@ protected: } }; + void parse_ECS_conductivity(rapidjson::Value &conductivity, CoolPropFluid & fluid) + { + fluid.transport.conductivity_ecs.reference_fluid = cpjson::get_string(conductivity,"reference_fluid"); + + // Parameters for correction polynomials + fluid.transport.conductivity_ecs.psi_a = cpjson::get_long_double_array(conductivity["psi"]["a"]); + fluid.transport.conductivity_ecs.psi_t = cpjson::get_long_double_array(conductivity["psi"]["t"]); + fluid.transport.conductivity_ecs.psi_rhomolar_reducing = cpjson::get_double(conductivity["psi"],"rhomolar_reducing"); + fluid.transport.conductivity_ecs.f_int_a = cpjson::get_long_double_array(conductivity["f_int"]["a"]); + fluid.transport.conductivity_ecs.f_int_t = cpjson::get_long_double_array(conductivity["f_int"]["t"]); + fluid.transport.conductivity_ecs.f_int_T_reducing = cpjson::get_double(conductivity["f_int"],"T_reducing"); + + fluid.transport.conductivity_using_ECS = true; + } + void parse_ECS_viscosity(rapidjson::Value &viscosity, CoolPropFluid & fluid) { fluid.transport.viscosity_ecs.reference_fluid = cpjson::get_string(viscosity,"reference_fluid"); @@ -446,7 +461,7 @@ protected: 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; + fluid.transport.viscosity_using_ECS = true; } /// Parse the transport properties @@ -635,6 +650,12 @@ protected: /// Parse the thermal conductivity data void parse_thermal_conductivity(rapidjson::Value &conductivity, CoolPropFluid & fluid) { + // If it is using ECS, set ECS parameters and quit + if (conductivity.HasMember("type") && !cpjson::get_string(conductivity, "type").compare("ECS")){ + parse_ECS_conductivity(conductivity, fluid); + return; + } + if (conductivity.HasMember("hardcoded")){ std::string target = cpjson::get_string(conductivity, "hardcoded"); if (!target.compare("Water")){ diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 4dcf839d..c444597f 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -199,9 +199,14 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void) CoolPropFluid &component = *(components[0]); // Check if using ECS - if (component.transport.using_ECS) + if (component.transport.viscosity_using_ECS) { - return TransportRoutines::viscosity_ECS(*this, *this); + // Get reference fluid name + std::string fluid_name = component.transport.viscosity_ecs.reference_fluid; + // Get a managed pointer to the reference fluid for ECS + std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(std::vector(1, fluid_name))); + // Get the viscosity using ECS + return TransportRoutines::viscosity_ECS(*this, *(ref_fluid.get())); } if (component.transport.hardcoded_viscosity != CoolProp::TransportPropertyData::VISCOSITY_NOT_HARDCODED) @@ -234,13 +239,42 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void) throw NotImplementedError(format("viscosity not implemented for mixtures")); } } +long double HelmholtzEOSMixtureBackend::calc_conductivity_background(void) +{ + // Residual part + long double lambda_residual = _HUGE; + switch(components[0]->transport.conductivity_residual.type) + { + case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL: + lambda_residual = TransportRoutines::conductivity_residual_polynomial(*this); break; + case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL_AND_EXPONENTIAL: + lambda_residual = TransportRoutines::conductivity_residual_polynomial_and_exponential(*this); break; + default: + throw ValueError(format("residual conductivity type [%d] is invalid for fluid %s", components[0]->transport.conductivity_residual.type, name().c_str())); + } + return lambda_residual; +} long double HelmholtzEOSMixtureBackend::calc_conductivity(void) { if (is_pure_or_pseudopure) { - if (components[0]->transport.hardcoded_conductivity != CoolProp::TransportPropertyData::CONDUCTIVITY_NOT_HARDCODED) + // Get a reference for code cleanness + CoolPropFluid &component = *(components[0]); + + // Check if using ECS + if (component.transport.conductivity_using_ECS) { - switch(components[0]->transport.hardcoded_conductivity) + // Get reference fluid name + std::string fluid_name = component.transport.conductivity_ecs.reference_fluid; + // Get a managed pointer to the reference fluid for ECS + std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(std::vector(1, fluid_name))); + // Get the viscosity using ECS + return TransportRoutines::conductivity_ECS(*this, *(ref_fluid.get())); + } + + if (component.transport.hardcoded_conductivity != CoolProp::TransportPropertyData::CONDUCTIVITY_NOT_HARDCODED) + { + switch(component.transport.hardcoded_conductivity) { case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_WATER: return TransportRoutines::conductivity_hardcoded_water(*this); @@ -255,7 +289,7 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void) // Dilute part long double lambda_dilute = _HUGE; - switch(components[0]->transport.conductivity_dilute.type) + switch(component.transport.conductivity_dilute.type) { case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_RATIO_POLYNOMIALS: lambda_dilute = TransportRoutines::conductivity_dilute_ratio_polynomials(*this); break; @@ -269,21 +303,11 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void) throw ValueError(format("dilute conductivity type [%d] is invalid for fluid %s", components[0]->transport.conductivity_dilute.type, name().c_str())); } - // Residual part - long double lambda_residual = _HUGE; - switch(components[0]->transport.conductivity_residual.type) - { - case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL: - lambda_residual = TransportRoutines::conductivity_residual_polynomial(*this); break; - case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL_AND_EXPONENTIAL: - lambda_residual = TransportRoutines::conductivity_residual_polynomial_and_exponential(*this); break; - default: - throw ValueError(format("residual conductivity type [%d] is invalid for fluid %s", components[0]->transport.conductivity_residual.type, name().c_str())); - } + long double lambda_residual = calc_conductivity_background(); // Critical part long double lambda_critical = _HUGE; - switch(components[0]->transport.conductivity_critical.type) + switch(component.transport.conductivity_critical.type) { case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_SIMPLIFIED_OLCHOWY_SENGERS: lambda_critical = TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(*this); break; @@ -1305,7 +1329,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double if (!ValidNumber(rhomolar)){throw ValueError();} return rhomolar; } - catch(std::exception &e) + catch(std::exception &) { try{ // Next we try with Secant method @@ -1313,7 +1337,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double if (!ValidNumber(rhomolar)){throw ValueError();} return rhomolar; } - catch(std::exception &e) + catch(std::exception &) { Secant(resid, rhomolar_guess, 0.0001*rhomolar_guess, 1e-8, 100, errstring); return _HUGE; @@ -1538,6 +1562,19 @@ long double HelmholtzEOSMixtureBackend::calc_cpmolar(void) return static_cast(_cpmolar); } +long double HelmholtzEOSMixtureBackend::calc_cpmolar_idealgas(void) +{ + // Calculate the reducing parameters + _delta = _rhomolar/_reducing.rhomolar; + _tau = _reducing.T/_T; + + // Calculate derivatives if needed, or just use cached values + long double d2a0_dTau2 = d2alpha0_dTau2(); + long double R_u = static_cast(_gas_constant); + + // Get cp of the ideal gas + return R_u*-pow(_tau.pt(),2)*d2a0_dTau2; +} long double HelmholtzEOSMixtureBackend::calc_speed_sound(void) { // Calculate the reducing parameters diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index d21339eb..8d133477 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -93,6 +93,7 @@ public: long double calc_pressure(void); long double calc_cvmolar(void); long double calc_cpmolar(void); + long double calc_cpmolar_idealgas(void); long double calc_hmolar(void); long double calc_smolar(void); long double calc_pressure_nocache(long double T, long double rhomolar); @@ -123,6 +124,7 @@ public: long double calc_viscosity_background(void); long double calc_viscosity_background(long double eta_dilute); long double calc_conductivity(void); + long double calc_conductivity_background(void); long double calc_Tmax(void); long double calc_pmax(void); diff --git a/src/Backends/Helmholtz/TransportRoutines.cpp b/src/Backends/Helmholtz/TransportRoutines.cpp index 3b76b6b5..07d24e1d 100644 --- a/src/Backends/Helmholtz/TransportRoutines.cpp +++ b/src/Backends/Helmholtz/TransportRoutines.cpp @@ -542,7 +542,7 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers( double Xref = Pcrit/pow(rhoc, 2)*HEOS.rhomolar()/dp_drho_ref*Tref/HEOS.T(); num = X - Xref; - // no critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff + // No critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff // See also Lemmon, IJT, 2004, page 27 if (num < DBL_EPSILON*10) return 0.0; @@ -827,21 +827,18 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H rhocmolar = HEOS.rhomolar_critical(), rhocmolar0 = HEOS_Reference.rhomolar_critical(); + // Get a reference to the ECS data CoolProp::ViscosityECSVariables &ECS = HEOS.components[0]->transport.viscosity_ecs; - - std::vector &a = ECS.psi_a, &t = ECS.psi_t; // 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.psi_rhomolar_reducing, t[i]); + for (std::size_t i=0; i < ECS.psi_a.size(); i++){ + psi += ECS.psi_a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, ECS.psi_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; @@ -869,4 +866,76 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H return eta; } + +long double TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference) +{ + // Collect some parameters + long double M = HEOS.molar_mass(), + M_kmol = M*1000, + M0 = HEOS_Reference.molar_mass(), + Tc = HEOS.T_critical(), + Tc0 = HEOS_Reference.T_critical(), + rhocmolar = HEOS.rhomolar_critical(), + rhocmolar0 = HEOS_Reference.rhomolar_critical(), + R_u = HEOS.gas_constant(), + R = HEOS.gas_constant()/HEOS.molar_mass(), //[J/kg/K] + R_kJkgK = R_u/M_kmol; + + // Get a reference to the ECS data + CoolProp::ConductivityECSVariables &ECS = HEOS.components[0]->transport.conductivity_ecs; + + // The correction polynomial psi_eta in rho/rho_red + double psi = 0; + for (std::size_t i=0; i < ECS.psi_a.size(); ++i){ + psi += ECS.psi_a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, ECS.psi_t[i]); + } + + // The correction polynomial f_int in T/T_red + double fint = 0; + for (std::size_t i=0; i < ECS.f_int_a.size(); ++i){ + fint += ECS.f_int_a[i]*pow(HEOS.T()/ECS.f_int_T_reducing, ECS.f_int_t[i]); + } + + // The dilute gas density for the fluid of interest [uPa-s] + long double eta_dilute = viscosity_dilute_kinetic_theory(HEOS)*1e6; + + // The mass specific ideal gas constant-pressure specific heat [J/kg/K] + long double cp0 = HEOS.calc_cpmolar_idealgas()/HEOS.molar_mass(); + + // The internal contribution to the thermal conductivity [W/m/K] + long double lambda_int = fint*eta_dilute*(cp0 - 2.5*R)/1e3; + + // The dilute gas contribution to the thermal conductivity [W/m/K] + long double lambda_dilute = 15.0e-3/4.0*R_kJkgK*eta_dilute; + + /// \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 conductivity [W/m/K] + long double lambda_resid = HEOS_Reference.calc_conductivity_background(); + + // The F factor + long double F_lambda = sqrt(f)*pow(h, static_cast(-2.0/3.0))*sqrt(M0/M); + + // The critical contribution from the fluid of interest [W/m/K] + long double lambda_critical = conductivity_critical_simplified_Olchowy_Sengers(HEOS); + + // The total thermal conductivity considering the contributions of the fluid of interest and the reference fluid [Pa-s] + long double lambda = lambda_int + lambda_dilute + lambda_resid*F_lambda + lambda_critical; + + return lambda; +} }; /* namespace CoolProp */ \ No newline at end of file diff --git a/src/Backends/Helmholtz/TransportRoutines.h b/src/Backends/Helmholtz/TransportRoutines.h index 3fe65fb9..ab298870 100644 --- a/src/Backends/Helmholtz/TransportRoutines.h +++ b/src/Backends/Helmholtz/TransportRoutines.h @@ -193,6 +193,8 @@ public: */ static long double viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference); + static long double conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference); + }; /* class TransportRoutines */ diff --git a/src/HumidAirProp.cpp b/src/HumidAirProp.cpp index 46d5fcd3..16242419 100644 --- a/src/HumidAirProp.cpp +++ b/src/HumidAirProp.cpp @@ -16,19 +16,18 @@ #include #include -CoolProp::AbstractStateWrapper Water; -CoolProp::AbstractStateWrapper Air; +std::tr1::shared_ptr Water, Air; namespace HumidAir { void check_fluid_instantiation() { - if (Water.empty()){ - Water.set("HEOS", "Water"); + if (!Water.get()){ + Water.reset(CoolProp::AbstractState::factory("HEOS", "Water")); } - if (Air.empty()){ - Air.set("HEOS", "Air"); + if (!Air.get()){ + Air.reset(CoolProp::AbstractState::factory("HEOS", "Air")); } }; @@ -42,60 +41,60 @@ double f_factor(double T, double p); static double MM_Air(void) { check_fluid_instantiation(); - return Air.keyed_output(CoolProp::imolar_mass); + return Air->keyed_output(CoolProp::imolar_mass); } static double MM_Water(void) { check_fluid_instantiation(); - return Water.keyed_output(CoolProp::imolar_mass); + return Water->keyed_output(CoolProp::imolar_mass); } static double B_Air(double T) { check_fluid_instantiation(); - Air.update(CoolProp::DmolarT_INPUTS,1e-12,T); - return Air.keyed_output(CoolProp::iBvirial); + Air->update(CoolProp::DmolarT_INPUTS,1e-12,T); + return Air->keyed_output(CoolProp::iBvirial); } static double dBdT_Air(double T) { check_fluid_instantiation(); - Air.update(CoolProp::DmolarT_INPUTS,1e-12,T); - return Air.keyed_output(CoolProp::idBvirial_dT); + Air->update(CoolProp::DmolarT_INPUTS,1e-12,T); + return Air->keyed_output(CoolProp::idBvirial_dT); } static double B_Water(double T) { check_fluid_instantiation(); - Water.update(CoolProp::DmolarT_INPUTS,1e-12,T); - return Water.keyed_output(CoolProp::iBvirial); + Water->update(CoolProp::DmolarT_INPUTS,1e-12,T); + return Water->keyed_output(CoolProp::iBvirial); } static double dBdT_Water(double T) { check_fluid_instantiation(); - Water.update(CoolProp::DmolarT_INPUTS,1e-12,T); - return Water.keyed_output(CoolProp::idBvirial_dT); + Water->update(CoolProp::DmolarT_INPUTS,1e-12,T); + return Water->keyed_output(CoolProp::idBvirial_dT); } static double C_Air(double T) { check_fluid_instantiation(); - Air.update(CoolProp::DmolarT_INPUTS,1e-12,T); - return Air.keyed_output(CoolProp::iCvirial); + Air->update(CoolProp::DmolarT_INPUTS,1e-12,T); + return Air->keyed_output(CoolProp::iCvirial); } static double dCdT_Air(double T) { check_fluid_instantiation(); - Air.update(CoolProp::DmolarT_INPUTS,1e-12,T); - return Air.keyed_output(CoolProp::idCvirial_dT); + Air->update(CoolProp::DmolarT_INPUTS,1e-12,T); + return Air->keyed_output(CoolProp::idCvirial_dT); } static double C_Water(double T) { check_fluid_instantiation(); - Water.update(CoolProp::DmolarT_INPUTS,1e-20,T); - return Water.keyed_output(CoolProp::iCvirial); + Water->update(CoolProp::DmolarT_INPUTS,1e-20,T); + return Water->keyed_output(CoolProp::iCvirial); } static double dCdT_Water(double T) { check_fluid_instantiation(); - Water.update(CoolProp::DmolarT_INPUTS,1e-12,T); - return Water.keyed_output(CoolProp::idCvirial_dT); + Water->update(CoolProp::DmolarT_INPUTS,1e-12,T); + return Water->keyed_output(CoolProp::idCvirial_dT); } void UseVirialCorrelations(int flag) { @@ -364,7 +363,7 @@ static double dC_m_dT(double T, double psi_w) // NDG for fluid EOS for virial terms Tj=132.6312; tau_Air=Tj/T; - tau_Water=Water.keyed_output(CoolProp::iT_reducing)/T; + tau_Water=Water->keyed_output(CoolProp::iT_reducing)/T; if (FlagUseVirialCorrelations) { dC_dT_aaa=-2.46582342273e-10 +4.425401935447e-12*T -3.669987371644e-14*pow(T,2) +1.765891183964e-16*pow(T,3) -5.240097805744e-19*pow(T,4) +9.502177003614e-22*pow(T,5) -9.694252610339e-25*pow(T,6) +4.276261986741e-28*pow(T,7); @@ -390,8 +389,8 @@ static double HenryConstant(double T) double p_ws,beta_N2,beta_O2,beta_Ar,beta_a,tau,Tr,Tc=647.096; Tr=T/Tc; tau=1-Tr; - Water.update(CoolProp::QT_INPUTS, 1.0, T); - p_ws = Water.keyed_output(CoolProp::iP); //[Pa] + Water->update(CoolProp::QT_INPUTS, 1.0, T); + p_ws = Water->keyed_output(CoolProp::iP); //[Pa] beta_N2=p_ws*exp(-9.67578/Tr+4.72162*pow(tau,0.355)/Tr+11.70585*pow(Tr,-0.41)*exp(tau)); beta_O2=p_ws*exp(-9.44833/Tr+4.43822*pow(tau,0.355)/Tr+11.42005*pow(Tr,-0.41)*exp(tau)); beta_Ar=p_ws*exp(-8.40954/Tr+4.29587*pow(tau,0.355)/Tr+10.52779*pow(Tr,-0.41)*exp(tau)); @@ -411,8 +410,8 @@ double isothermal_compressibility(double T, double p) } else { - Water.update(CoolProp::PT_INPUTS, p, T); - k_T = Water.keyed_output(CoolProp::iisothermal_compressibility); + Water->update(CoolProp::PT_INPUTS, p, T); + k_T = Water->keyed_output(CoolProp::iisothermal_compressibility); } } else @@ -436,8 +435,8 @@ double f_factor(double T, double p) // It is liquid water p_ws=CoolProp::PropsSI("P","T",T,"Q",0,"Water"); beta_H = HenryConstant(T); //[1/Pa] - Water.update(CoolProp::PT_INPUTS, p, T); - vbar_ws = 1.0/Water.keyed_output(CoolProp::iDmolar); //[m^3/mol] + Water->update(CoolProp::PT_INPUTS, p, T); + vbar_ws = 1.0/Water->keyed_output(CoolProp::iDmolar); //[m^3/mol] } else { @@ -459,7 +458,7 @@ double f_factor(double T, double p) // NDG for fluid EOS for virial terms Tj=132.6312; tau_Air=Tj/T; - tau_Water=Water.keyed_output(CoolProp::iT_reducing)/T; + tau_Water=Water->keyed_output(CoolProp::iT_reducing)/T; if (FlagUseVirialCorrelations) { B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2) @@ -560,11 +559,11 @@ double Viscosity(double T, double p, double psi_w) Mw=MM_Water(); Ma=MM_Air(); // Viscosity of dry air at dry-bulb temp and total pressure - Air.update(CoolProp::PT_INPUTS,p,T); - mu_a=Air.keyed_output(CoolProp::iviscosity); + Air->update(CoolProp::PT_INPUTS,p,T); + mu_a=Air->keyed_output(CoolProp::iviscosity); // Viscosity of pure saturated water at dry-bulb temperature - Water.update(CoolProp::PQ_INPUTS,p,1); - mu_w=Water.keyed_output(CoolProp::iviscosity); + Water->update(CoolProp::PQ_INPUTS,p,1); + mu_w=Water->keyed_output(CoolProp::iviscosity); Phi_av=sqrt(2.0)/4.0*pow(1+Ma/Mw,-0.5)*pow(1+sqrt(mu_a/mu_w)*pow(Mw/Ma,0.25),2); //[-] Phi_va=sqrt(2.0)/4.0*pow(1+Mw/Ma,-0.5)*pow(1+sqrt(mu_w/mu_a)*pow(Ma/Mw,0.25),2); //[-] return (1-psi_w)*mu_a/((1-psi_w)+psi_w*Phi_av)+psi_w*mu_w/(psi_w+(1-psi_w)*Phi_va); @@ -583,13 +582,13 @@ double Conductivity(double T, double p, double psi_w) Ma=MM_Air(); // Viscosity of dry air at dry-bulb temp and total pressure - Air.update(CoolProp::PT_INPUTS,p,T); - mu_a=Air.keyed_output(CoolProp::iviscosity); - k_a=Air.keyed_output(CoolProp::iconductivity); + Air->update(CoolProp::PT_INPUTS,p,T); + mu_a=Air->keyed_output(CoolProp::iviscosity); + k_a=Air->keyed_output(CoolProp::iconductivity); // Viscosity of pure saturated water at dry-bulb temperature - Water.update(CoolProp::PQ_INPUTS,p,1); - mu_w=Water.keyed_output(CoolProp::iviscosity); - k_w=Water.keyed_output(CoolProp::iconductivity); + Water->update(CoolProp::PQ_INPUTS,p,1); + mu_w=Water->keyed_output(CoolProp::iviscosity); + k_w=Water->keyed_output(CoolProp::iconductivity); Phi_av=sqrt(2.0)/4.0*pow(1+Ma/Mw,-0.5)*pow(1+sqrt(mu_a/mu_w)*pow(Mw/Ma,0.25),2); //[-] Phi_va=sqrt(2.0)/4.0*pow(1+Mw/Ma,-0.5)*pow(1+sqrt(mu_w/mu_a)*pow(Ma/Mw,0.25),2); //[-] return (1-psi_w)*k_a/((1-psi_w)+psi_w*Phi_av)+psi_w*k_w/(psi_w+(1-psi_w)*Phi_va); @@ -637,19 +636,19 @@ double IdealGasMolarEnthalpy_Water(double T, double vmolar) double hbar_w_0, tau, rhomolar, hbar_w; // Ideal-Gas contribution to enthalpy of water hbar_w_0 = -0.01102303806;//[J/mol] - tau = Water.keyed_output(CoolProp::iT_reducing)/T; + tau = Water->keyed_output(CoolProp::iT_reducing)/T; rhomolar = 1/vmolar; //[mol/m^3] - Water.update(CoolProp::DmolarT_INPUTS, rhomolar, T); - hbar_w = hbar_w_0+R_bar*T*(1+tau*Water.keyed_output(CoolProp::idalpha0_dtau_constdelta)); + Water->update(CoolProp::DmolarT_INPUTS, rhomolar, T); + hbar_w = hbar_w_0+R_bar*T*(1+tau*Water->keyed_output(CoolProp::idalpha0_dtau_constdelta)); return hbar_w; } double IdealGasMolarEntropy_Water(double T, double p) { double sbar_w, tau, R_bar; R_bar = 8.314371; //[J/mol/K] - tau = Water.keyed_output(CoolProp::iT_reducing)/T; - Water.update(CoolProp::DmolarT_INPUTS,p/(R_bar*T),T); - sbar_w = R_bar*(tau*Water.keyed_output(CoolProp::idalpha0_dtau_constdelta)-Water.keyed_output(CoolProp::ialpha0)); //[kJ/kmol/K] + tau = Water->keyed_output(CoolProp::iT_reducing)/T; + Water->update(CoolProp::DmolarT_INPUTS,p/(R_bar*T),T); + sbar_w = R_bar*(tau*Water->keyed_output(CoolProp::idalpha0_dtau_constdelta)-Water->keyed_output(CoolProp::ialpha0)); //[kJ/kmol/K] return sbar_w; } double IdealGasMolarEnthalpy_Air(double T, double vmolar) @@ -661,9 +660,9 @@ double IdealGasMolarEnthalpy_Air(double T, double vmolar) tau = 132.6312/T; rhomolar = 1/vmolar; //[mol/m^3] R_bar_Lemmon = 8.314510; //[J/mol/K] - Air.update(CoolProp::DmolarT_INPUTS, rhomolar, T); - double dd = Air.keyed_output(CoolProp::idalpha0_dtau_constdelta); - hbar_a = hbar_a_0 + R_bar_Lemmon*T*(1+tau*Air.keyed_output(CoolProp::idalpha0_dtau_constdelta)); //[J/mol] + Air->update(CoolProp::DmolarT_INPUTS, rhomolar, T); + double dd = Air->keyed_output(CoolProp::idalpha0_dtau_constdelta); + hbar_a = hbar_a_0 + R_bar_Lemmon*T*(1+tau*Air->keyed_output(CoolProp::idalpha0_dtau_constdelta)); //[J/mol] return hbar_a; } double IdealGasMolarEntropy_Air(double T, double vmolar_a) @@ -678,9 +677,9 @@ double IdealGasMolarEntropy_Air(double T, double vmolar_a) vmolar_a_0 = R_bar_Lemmon*T0/p0; //[m^3/mol] - Air.update(CoolProp::DmolarT_INPUTS,1/vmolar_a_0,T); + Air->update(CoolProp::DmolarT_INPUTS,1/vmolar_a_0,T); - sbar_a=sbar_0_Lem+R_bar_Lemmon*(tau*Air.keyed_output(CoolProp::idalpha0_dtau_constdelta)-Air.keyed_output(CoolProp::ialpha0))+R_bar_Lemmon*log(vmolar_a/vmolar_a_0); //[J/mol/K] + sbar_a=sbar_0_Lem+R_bar_Lemmon*(tau*Air->keyed_output(CoolProp::idalpha0_dtau_constdelta)-Air->keyed_output(CoolProp::ialpha0))+R_bar_Lemmon*log(vmolar_a/vmolar_a_0); //[J/mol/K] return sbar_a; //[J/mol/K] } @@ -802,8 +801,8 @@ double DewpointTemperature(double T, double p, double psi_w) // 0.61165... is the triple point pressure of water in kPa if (p_w > 0.6116547241637944){ - Water.update(CoolProp::PQ_INPUTS, p, 1.0); - T0 = Water.keyed_output(CoolProp::iT); + Water->update(CoolProp::PQ_INPUTS, p, 1.0); + T0 = Water->keyed_output(CoolProp::iT); } else{ T0 = 268; @@ -821,8 +820,8 @@ double DewpointTemperature(double T, double p, double psi_w) if (Tdp >= 273.16) { // Saturation pressure at dewpoint [kPa] - Water.update(CoolProp::QT_INPUTS, 0.0, Tdp); - p_ws_dp = Water.keyed_output(CoolProp::iP); + Water->update(CoolProp::QT_INPUTS, 0.0, Tdp); + p_ws_dp = Water->keyed_output(CoolProp::iP); } else { @@ -872,8 +871,8 @@ public: if (Twb > 273.16) { // Saturation pressure at wetbulb temperature [Pa] - Water.update(CoolProp::QT_INPUTS,0,Twb); - p_ws_wb= Water.keyed_output(CoolProp::iP); + Water->update(CoolProp::QT_INPUTS,0,Twb); + p_ws_wb= Water->keyed_output(CoolProp::iP); } else { @@ -890,8 +889,8 @@ public: if (Twb > 273.16) { // Enthalpy of water [J/kg_water] - Water.update(CoolProp::PT_INPUTS, _p, Twb); - h_w = Water.keyed_output(CoolProp::iHmass); //[J/kg_water] + Water->update(CoolProp::PT_INPUTS, _p, Twb); + h_w = Water->keyed_output(CoolProp::iHmass); //[J/kg_water] } else { @@ -939,8 +938,8 @@ double WetbulbTemperature(double T, double p, double psi_w) // If the temperature is above the saturation temperature corresponding to the atmospheric pressure, // then the maximum value for the wetbulb temperature is the saturation temperature double Tmax = T; - Water.update(CoolProp::PQ_INPUTS,p,1.0); - double Tsat = Water.keyed_output(CoolProp::iT); + Water->update(CoolProp::PQ_INPUTS,p,1.0); + double Tsat = Water->keyed_output(CoolProp::iT); if (T >= Tsat) { Tmax = Tsat; @@ -1031,8 +1030,8 @@ double MoleFractionWater(double T, double p, int HumInput, double InVal) if (T>=273.16) { // Saturation pressure [Pa] - Water.update(CoolProp::QT_INPUTS,0,T); - p_ws= Water.keyed_output(CoolProp::iP);; + Water->update(CoolProp::QT_INPUTS,0,T); + p_ws= Water->keyed_output(CoolProp::iP);; } else { @@ -1055,8 +1054,8 @@ double MoleFractionWater(double T, double p, int HumInput, double InVal) // Saturation pressure at dewpoint [Pa] if (Tdp>=273.16) { - Water.update(CoolProp::QT_INPUTS,0,Tdp); - p_ws_dp = Water.keyed_output(CoolProp::iP); //[Pa] + Water->update(CoolProp::QT_INPUTS,0,Tdp); + p_ws_dp = Water->keyed_output(CoolProp::iP); //[Pa] } else{ // Sublimation pressure [Pa] @@ -1081,8 +1080,8 @@ double RelativeHumidity(double T, double p, double psi_w) double p_ws, f, p_s, W; if (T >= 273.16){ // Saturation pressure [Pa] - Water.update(CoolProp::QT_INPUTS, 0, T); - p_ws = Water.keyed_output(CoolProp::iP); //[Pa] + Water->update(CoolProp::QT_INPUTS, 0, T); + p_ws = Water->keyed_output(CoolProp::iP); //[Pa] } else{ // sublimation pressure [Pa] @@ -1400,7 +1399,7 @@ double HAProps_Aux(const char* Name,double T, double p, double W, char *units) Tj=132.6312; tau_Air=Tj/T; - tau_Water=Water.keyed_output(CoolProp::iT_critical)/T; + tau_Water=Water->keyed_output(CoolProp::iT_critical)/T; try{ if (!strcmp(Name,"Baa")) @@ -1497,8 +1496,8 @@ double HAProps_Aux(const char* Name,double T, double p, double W, char *units) strcpy(units,"1/Pa"); if (T>273.16) { - Water.update(CoolProp::PT_INPUTS, p, T); - return Water.keyed_output(CoolProp::iisothermal_compressibility); + Water->update(CoolProp::PT_INPUTS, p, T); + return Water->keyed_output(CoolProp::iisothermal_compressibility); } else return IsothermCompress_Ice(T,p); //[1/Pa] @@ -1508,8 +1507,8 @@ double HAProps_Aux(const char* Name,double T, double p, double W, char *units) strcpy(units,"Pa"); if (T>273.16) { - Water.update(CoolProp::QT_INPUTS, 0, T); - return Water.keyed_output(CoolProp::iP); + Water->update(CoolProp::QT_INPUTS, 0, T); + return Water->keyed_output(CoolProp::iP); } else return psub_Ice(T); @@ -1519,8 +1518,8 @@ double HAProps_Aux(const char* Name,double T, double p, double W, char *units) strcpy(units,"m^3/mol"); if (T>273.16) { - Water.update(CoolProp::QT_INPUTS, 0, T); - return 1.0/Water.keyed_output(CoolProp::iDmolar); + Water->update(CoolProp::QT_INPUTS, 0, T); + return 1.0/Water->keyed_output(CoolProp::iDmolar); } else {