From edf39c7387b9c8489fbced06397d109c07c830e5 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Mon, 8 Sep 2014 17:34:05 +0200 Subject: [PATCH] Finished normalization of the gas constants Signed-off-by: Ian Bell --- include/CoolPropFluid.h | 3 +- src/Backends/Helmholtz/Configuration.cpp | 42 +++++++ src/Backends/Helmholtz/Configuration.h | 27 +++++ src/Backends/Helmholtz/Fluids/FluidLibrary.h | 108 +++++++++++++++--- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 28 +++-- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 1 + src/Backends/Helmholtz/MixtureDerivatives.h | 4 +- src/Backends/Helmholtz/ReducingFunctions.h | 18 ++- src/Backends/Helmholtz/TransportRoutines.cpp | 25 ++-- src/Backends/Helmholtz/VLERoutines.cpp | 11 +- 10 files changed, 224 insertions(+), 43 deletions(-) create mode 100644 src/Backends/Helmholtz/Configuration.cpp create mode 100644 src/Backends/Helmholtz/Configuration.h diff --git a/include/CoolPropFluid.h b/include/CoolPropFluid.h index eca84b00..8760f2f4 100644 --- a/include/CoolPropFluid.h +++ b/include/CoolPropFluid.h @@ -320,7 +320,8 @@ public: max_sat_T, ///< The state at the maximum saturation temperature for pseudo-pure max_sat_p; ///< The state at the maximum saturation pressure for pseudo-pure EOSLimits limits; ///< Limits on the EOS - double R_u, ///< The universal gas constant used for this EOS (usually, but not always, 8.314472 J/mol/K) + double R_u_specified, ///< The universal gas constant that was specified in the equation of state + R_u, ///< The universal gas constant used for this EOS (usually, but not always, 8.314472 J/mol/K) molar_mass, ///< The molar mass in kg/mol (note NOT kg/kmol) accentric, ///< The accentric factor \f$ \omega = -log_{10}\left(\frac{p_s(T/T_c=0.7)}{p_c}\right)-1\f$ Ttriple, ///< Triple point temperature (K) diff --git a/src/Backends/Helmholtz/Configuration.cpp b/src/Backends/Helmholtz/Configuration.cpp new file mode 100644 index 00000000..0c209124 --- /dev/null +++ b/src/Backends/Helmholtz/Configuration.cpp @@ -0,0 +1,42 @@ +#include "Configuration.h" + +namespace CoolProp +{ +// +//Configuration::Configuration() +//{ +//} +// +//Configuration::~Configuration() +//{ +//} + +bool get_config_bool(configuration_keys key) +{ + switch(key) + { + case NORMALIZE_GAS_CONSTANTS: + return true; + default: + throw ValueError(format("%d is invalid key to get_config_bool",key)); + } +} +double get_config_double(configuration_keys key) +{ + switch(key) + { + default: + throw ValueError(format("%d is invalid key to get_config_double",key)); + } +} +std::string get_config_string(configuration_keys key) +{ + switch(key) + { + default: + throw ValueError(format("%d is invalid key to get_config_string",key)); + } +} + +} + diff --git a/src/Backends/Helmholtz/Configuration.h b/src/Backends/Helmholtz/Configuration.h new file mode 100644 index 00000000..b061ebde --- /dev/null +++ b/src/Backends/Helmholtz/Configuration.h @@ -0,0 +1,27 @@ +#ifndef COOLPROP_CONFIGURATION +#define COOLPROP_CONFIGURATION + +#include "Exceptions.h" +#include "CoolPropTools.h" + +enum configuration_keys {NORMALIZE_GAS_CONSTANTS}; + +namespace CoolProp +{ +// +//class Configuration +//{ +//public: +// Configuration(); +// ~Configuration(); +// +//}; + +/// Return the value of a boolean key from the configuration +bool get_config_bool(configuration_keys key); +double get_config_double(configuration_keys key); +std::string get_config_string(configuration_keys key); + +} + +#endif // COOLPROP_CONFIGURATION diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index b4d0cbe8..cbeafe80 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -8,6 +8,7 @@ #include #include +#include "../Configuration.h" namespace CoolProp{ @@ -49,8 +50,10 @@ protected: assert(n.size() == d.size()); assert(n.size() == t.size()); assert(n.size() == l.size()); - // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures - for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } EOS.alphar.GenExp.add_Power(n,d,t,l); } else if (!type.compare("ResidualHelmholtzGaussian")) @@ -69,8 +72,10 @@ protected: assert(n.size() == beta.size()); assert(n.size() == gamma.size()); - // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures - for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } EOS.alphar.GenExp.add_Gaussian(n,d,t,eta,epsilon,beta,gamma); } @@ -93,8 +98,10 @@ protected: assert(n.size() == C.size()); assert(n.size() == D.size()); - // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures - for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } EOS.alphar.NonAnalytic = ResidualHelmholtzNonAnalytic(n,a,b,beta,A,B,C,D); } @@ -110,8 +117,10 @@ protected: assert(n.size() == l.size()); assert(n.size() == m.size()); - // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures - for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } EOS.alphar.GenExp.add_Lemmon2005(n,d,t,l,m); } @@ -127,8 +136,10 @@ protected: assert(n.size() == g.size()); assert(n.size() == l.size()); - // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures - for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } EOS.alphar.GenExp.add_Exponential(n,d,t,g,l); } @@ -140,8 +151,11 @@ protected: long double epsilonbar = cpjson::get_double(contribution,"epsilonbar"); long double vbarn = cpjson::get_double(contribution,"vbarn"); long double kappabar = cpjson::get_double(contribution,"kappabar"); - // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures - a *= EOS.R_u/R_u_CODATA; + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + a *= EOS.R_u/R_u_CODATA; + } EOS.alphar.SAFT = ResidualHelmholtzSAFTAssociating(a,m,epsilonbar,vbarn,kappabar); } else @@ -171,6 +185,13 @@ protected: if (EOS.alpha0.Lead.is_enabled() == true){throw ValueError("Cannot add ");} long double a1 = cpjson::get_double(contribution,"a1"); long double a2 = cpjson::get_double(contribution,"a2"); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + a1 *= EOS.R_u/R_u_CODATA; + a2 *= EOS.R_u/R_u_CODATA; + } + EOS.alpha0.Lead = IdealHelmholtzLead(a1, a2); } else if (!type.compare("IdealGasHelmholtzPower")) @@ -178,12 +199,24 @@ protected: if (EOS.alpha0.Power.is_enabled() == true){throw ValueError("Cannot add ");} std::vector n = cpjson::get_long_double_array(contribution["n"]); std::vector t = cpjson::get_long_double_array(contribution["t"]); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } + EOS.alpha0.Power = IdealHelmholtzPower(n, t); } else if (!type.compare("IdealGasHelmholtzLogTau")) { if (EOS.alpha0.LogTau.is_enabled() == true){throw ValueError("Cannot add ");} long double a = cpjson::get_double(contribution,"a"); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + a *= EOS.R_u/R_u_CODATA; + } + EOS.alpha0.LogTau = IdealHelmholtzLogTau(a); } else if (!type.compare("IdealGasHelmholtzPlanckEinsteinGeneralized")) @@ -194,6 +227,12 @@ protected: std::vector c = cpjson::get_long_double_array(contribution["c"]); std::vector d = cpjson::get_long_double_array(contribution["d"]); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } + if (EOS.alpha0.PlanckEinstein.is_enabled() == true){ EOS.alpha0.PlanckEinstein.extend(n, t, c, d); } @@ -210,6 +249,12 @@ protected: for (std::size_t i = 0; i < t.size(); ++i){ t[i] *= -1;} std::vector c(n.size(), 1); std::vector d(c.size(), -1); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } + if (EOS.alpha0.PlanckEinstein.is_enabled() == true){ EOS.alpha0.PlanckEinstein.extend(n, t, c, d); } @@ -223,6 +268,11 @@ protected: long double cp_over_R = cpjson::get_double(contribution, "cp_over_R"); long double Tc = cpjson::get_double(contribution, "Tc"); long double T0 = cpjson::get_double(contribution, "T0"); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + cp_over_R *= EOS.R_u/R_u_CODATA; + } + EOS.alpha0.CP0Constant = IdealHelmholtzCP0Constant(cp_over_R, Tc, T0); } else if (!type.compare("IdealGasHelmholtzCP0PolyT")) @@ -232,6 +282,12 @@ protected: std::vector t = cpjson::get_long_double_array(contribution["t"]); long double Tc = cpjson::get_double(contribution, "Tc"); long double T0 = cpjson::get_double(contribution, "T0"); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < c.size(); ++i){ c[i] *= EOS.R_u/R_u_CODATA; } + } + EOS.alpha0.CP0PolyT = IdealHelmholtzCP0PolyT(c, t, Tc, T0); } else if (!type.compare("IdealGasHelmholtzCP0AlyLee")) @@ -244,6 +300,11 @@ protected: // Take the constant term if nonzero and set it as a polyT term if (std::abs(constants[0]) > 1e-14){ std::vector c(1,constants[0]), t(1,0); + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < c.size(); ++i){ c[i] *= EOS.R_u/R_u_CODATA; } + } + if (EOS.alpha0.CP0PolyT.is_enabled() == true){ EOS.alpha0.CP0PolyT.extend(c,t); } @@ -268,6 +329,11 @@ protected: c.push_back(1); d.push_back(1); } + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + for (std::size_t i = 0; i < n.size(); ++i){ n[i] *= EOS.R_u/R_u_CODATA; } + } if (EOS.alpha0.PlanckEinstein.is_enabled() == true){ EOS.alpha0.PlanckEinstein.extend(n, t, c, d); @@ -281,6 +347,13 @@ protected: long double a1 = cpjson::get_double(contribution, "a1"); long double a2 = cpjson::get_double(contribution, "a2"); std::string reference = cpjson::get_string(contribution, "reference"); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Normalize the gas constants to yield the same (CODATA 2010) value for all fluids - needed for mixtures + a1 *= EOS.R_u/R_u_CODATA; + a2 *= EOS.R_u/R_u_CODATA; + } + EOS.alpha0.EnthalpyEntropyOffset = IdealHelmholtzEnthalpyEntropyOffset(a1, a2, reference); } else @@ -386,10 +459,13 @@ protected: // Validate the equation of state that was just created EOS.validate(); - // Set the universal gas constant to the CODATA 2010 value for consistency. - // This must be the LAST step since the loaded value is used to adjust the coefficients n_i - EOS.R_u = R_u_CODATA; - + // Set the specified gas constant value - this is saved separately so that if normalization is applied, this term can be used in pressure and enthalpy terms + EOS.R_u_specified = EOS.R_u; + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + // Set the universal gas constant to the CODATA 2010 value for consistency. + // This must be the LAST step since the loaded value is used to adjust the coefficients + EOS.R_u = R_u_CODATA; + } } /// Parse the list of possible equations of state diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 8db1965c..ea255909 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -172,6 +172,15 @@ long double HelmholtzEOSMixtureBackend::calc_gas_constant(void) } return summer; } +long double HelmholtzEOSMixtureBackend::calc_gas_constant_specified(void) +{ + double summer = 0; + for (unsigned int i = 0; i < components.size(); ++i) + { + summer += mole_fractions[i]*components[i]->pEOS->R_u_specified; + } + return summer; +} long double HelmholtzEOSMixtureBackend::calc_molar_mass(void) { double summer = 0; @@ -1288,7 +1297,8 @@ void get_dT_drho(HelmholtzEOSMixtureBackend *HEOS, parameters index, long double dT_dtau = -pow(T, 2)/Tr, R = HEOS->gas_constant(), delta = rho/rhor, - tau = Tr/T; + tau = Tr/T, + R_u_ratio = HEOS->calc_gas_constant_specified()/HEOS->gas_constant(); // This term falls out as a result of normalizing gas constants - it is R_u given by EOS divided by CODATA value; switch (index) { @@ -1299,9 +1309,9 @@ void get_dT_drho(HelmholtzEOSMixtureBackend *HEOS, parameters index, long double case iP: { // dp/drho|T - drho = R*T*(1+2*delta*HEOS->dalphar_dDelta()+pow(delta, 2)*HEOS->d2alphar_dDelta2()); + drho = R*T*(R_u_ratio+2*delta*HEOS->dalphar_dDelta()+pow(delta, 2)*HEOS->d2alphar_dDelta2()); // dp/dT|rho - dT = rho*R*(1+delta*HEOS->dalphar_dDelta() - tau*delta*HEOS->d2alphar_dDelta_dTau()); + dT = rho*R*(R_u_ratio+delta*HEOS->dalphar_dDelta() - tau*delta*HEOS->d2alphar_dDelta_dTau()); break; } case iHmolar: @@ -1789,9 +1799,10 @@ long double HelmholtzEOSMixtureBackend::calc_pressure(void) // Calculate derivative if needed long double dar_dDelta = dalphar_dDelta(); long double R_u = gas_constant(); + long double R_u_ratio = calc_gas_constant_specified()/R_u; // This term falls out as a result of normalizing gas constants - it is R_u given by EOS divided by CODATA value // Get pressure - _p = _rhomolar*R_u*_T*(1+_delta.pt()*dar_dDelta); + _p = _rhomolar*R_u*_T*(R_u_ratio + _delta.pt()*dar_dDelta); //std::cout << format("p: %13.12f %13.12f %10.9f %10.9f %10.9f %10.9f %g\n",_T,_rhomolar,_tau,_delta,mole_fractions[0],dar_dDelta,_p); //if (_p < 0){ @@ -1834,9 +1845,10 @@ long double HelmholtzEOSMixtureBackend::calc_hmolar(void) long double dar_dTau = dalphar_dTau(); long double dar_dDelta = dalphar_dDelta(); long double R_u = gas_constant(); + long double R_u_ratio = calc_gas_constant_specified()/R_u; // This term falls out as a result of normalizing gas constants - it is R_u given by EOS divided by CODATA value // Get molar enthalpy - _hmolar = R_u*_T*(1 + _tau.pt()*(da0_dTau+dar_dTau) + _delta.pt()*dar_dDelta); + _hmolar = R_u*_T*(R_u_ratio + _tau.pt()*(da0_dTau+dar_dTau) + _delta.pt()*dar_dDelta); return static_cast(_hmolar); } @@ -1960,9 +1972,10 @@ long double HelmholtzEOSMixtureBackend::calc_cpmolar(void) long double d2ar_dDelta_dTau = d2alphar_dDelta_dTau(); long double d2ar_dTau2 = d2alphar_dTau2(); long double R_u = static_cast(_gas_constant); + long double R_u_ratio = calc_gas_constant_specified()/R_u; // This term falls out as a result of normalizing gas constants - it is R_u given by EOS divided by CODATA value // Get cp - _cpmolar = R_u*(-pow(_tau.pt(),2)*(d2ar_dTau2 + d2a0_dTau2)+pow(1+_delta.pt()*dar_dDelta-_delta.pt()*_tau.pt()*d2ar_dDelta_dTau,2)/(1+2*_delta.pt()*dar_dDelta+pow(_delta.pt(),2)*d2ar_dDelta2)); + _cpmolar = R_u*(-pow(_tau.pt(),2)*(d2ar_dTau2 + d2a0_dTau2)+pow(R_u_ratio+_delta.pt()*dar_dDelta-_delta.pt()*_tau.pt()*d2ar_dDelta_dTau,2)/(R_u_ratio+2*_delta.pt()*dar_dDelta+pow(_delta.pt(),2)*d2ar_dDelta2)); return static_cast(_cpmolar); } @@ -2017,9 +2030,10 @@ long double HelmholtzEOSMixtureBackend::calc_gibbsmolar(void) long double a0 = alpha0(); long double dar_dDelta = dalphar_dDelta(); long double R_u = gas_constant(); + long double R_u_ratio = calc_gas_constant_specified()/R_u; // This term falls out as a result of normalizing gas constants - it is R_u given by EOS divided by CODATA value // Get molar gibbs function - _gibbsmolar = R_u*_T*(1 + a0 + ar +_delta.pt()*dar_dDelta); + _gibbsmolar = R_u*_T*(R_u_ratio + a0 + ar +_delta.pt()*dar_dDelta); return static_cast(_gibbsmolar); } diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 7b13c73f..c6cd0f25 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -112,6 +112,7 @@ public: long double calc_molar_mass(void); long double calc_gas_constant(void); + long double calc_gas_constant_specified(void); long double calc_Bvirial(void); long double calc_Cvirial(void); diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index 12a8022c..521357c7 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -97,7 +97,7 @@ class MixtureDerivatives{ * * From Witzke, Eqn. 3.14 * \f[ - * \left(\frac{\partial \ln(f_i)}{\partial T} \right)_{\rho,x} = -\frac{1}{T}\left(1-\tau\alphar^r_{\tau}-\tau n\left(\frac{\partial\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}}{\partial \tau}\right)_{\delta,\bar x} \right) + * \left(\frac{\partial \ln(f_i)}{\partial T} \right)_{\rho,x} = -\frac{1}{T}\left(1-\tau\alpha^r_{\tau}-\tau n\left(\frac{\partial\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}}{\partial \tau}\right)_{\delta,\bar x} \right) * \f] * @param HEOS The HelmholtzEOSMixtureBackend to be used */ @@ -107,7 +107,7 @@ class MixtureDerivatives{ * * From Witzke, Eqn. 3.15 * \f[ - * \left(\frac{\partial \ln(f_i)}{\partial \rho} \right)_{T, x} = \frac{1}{\rho}\left(1+\delta\alphar^r_{\delta}+\delta n\left(\frac{\partial\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}}{\partial \delta}\right)_{\tau,\bar x} \right) + * \left(\frac{\partial \ln(f_i)}{\partial \rho} \right)_{T, x} = \frac{1}{\rho}\left(1+\delta\alpha^r_{\delta}+\delta n\left(\frac{\partial\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}}{\partial \delta}\right)_{\tau,\bar x} \right) * \f] * @param HEOS The HelmholtzEOSMixtureBackend to be used */ diff --git a/src/Backends/Helmholtz/ReducingFunctions.h b/src/Backends/Helmholtz/ReducingFunctions.h index cc19888b..7c8781ab 100644 --- a/src/Backends/Helmholtz/ReducingFunctions.h +++ b/src/Backends/Helmholtz/ReducingFunctions.h @@ -53,16 +53,28 @@ public: /** \brief GERG 2004 Monograph equation 7.56: * + * If the \f$x_i\f$ are all independent * \f[ - * \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=1}^Nx_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right) + * \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right) + * \f] + * If \f$x_N = 1-\sum x_i\f$: + * \f[ + * \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right) * \f] */ long double d_ndTrdni_dxj__constxi(const std::vector &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); - /** \brief GERG 2004 Monograph equation 7.55: + /** \brief * + * GERG 2004 Monograph equation 7.55: + * If the \f$x_i\f$ are all independent * \f[ - * \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=1}^Nx_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right) + * \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right) * \f] + * Gernert, JPCRD, 2014, A28 + * If \f$x_N = 1-\sum x_i\f$: + * \f[ + * \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-2}x_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right) + * \f] */ long double d_ndrhorbardni_dxj__constxi(const std::vector &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); diff --git a/src/Backends/Helmholtz/TransportRoutines.cpp b/src/Backends/Helmholtz/TransportRoutines.cpp index a8b780bb..fcb492a9 100644 --- a/src/Backends/Helmholtz/TransportRoutines.cpp +++ b/src/Backends/Helmholtz/TransportRoutines.cpp @@ -247,7 +247,7 @@ long double TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBack { double x_mu=0.068,qc=1/1.9,qd=1/1.1,nu=0.630,gamma=1.239,zeta_0=0.13,LAMBDA_0=0.06,Tbar_R=1.5, pstar, Tstar, rhostar; double delta,tau,mubar_0,mubar_1,mubar_2,drhodp,drhodp_R,DeltaChibar,zeta,w,L,Y,psi_D,Tbar,rhobar; - double drhobar_dpbar,drhobar_dpbar_R,R_Water; + double drhobar_dpbar,drhobar_dpbar_R,R_Water,R_spec,R_ratio; pstar = 22.064e6; // [Pa] Tstar = 647.096; // [K] @@ -255,6 +255,8 @@ long double TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBack Tbar = HEOS.T()/Tstar; rhobar = HEOS.rhomass()/rhostar; R_Water = HEOS.gas_constant()/HEOS.molar_mass(); // [J/kg/K] + R_spec = HEOS.calc_gas_constant_specified()/HEOS.molar_mass(); //[J/kg/K] + R_ratio = R_spec/R_Water; // Dilute and finite gas portions visc_Helper(Tbar, rhobar, &mubar_0, &mubar_1); @@ -265,11 +267,11 @@ long double TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBack delta=rhobar; // "Normal" calculation tau=1/Tbar; - drhodp=1/(R_Water*HEOS.T()*(1+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2())); + drhodp=1/(R_Water*HEOS.T()*(R_ratio+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2())); drhobar_dpbar = pstar/rhostar*drhodp; // "Reducing" calculation tau=1/Tbar_R; - drhodp_R=1/(R_Water*Tbar_R*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau, delta))); + drhodp_R=1/(R_Water*Tbar_R*Tstar*(R_ratio+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau, delta))); drhobar_dpbar_R = pstar/rhostar*drhodp_R; DeltaChibar=rhobar*(drhobar_dpbar-drhobar_dpbar_R*Tbar_R/Tbar); @@ -627,7 +629,8 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers( Pcrit = HEOS.get_reducing_state().p, // [Pa] Tref, // [K] cp,cv,delta,num,zeta,mu,pi=M_PI, - OMEGA_tilde,OMEGA_tilde0; + OMEGA_tilde,OMEGA_tilde0, + R_ratio = HEOS.calc_gas_constant_specified()/HEOS.gas_constant(); if (ValidNumber(data.T_ref)) Tref = data.T_ref; @@ -636,11 +639,11 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers( delta = HEOS.delta(); - double dp_drho = HEOS.gas_constant()*HEOS.T()*(1+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2()); + double dp_drho = HEOS.gas_constant()*HEOS.T()*(R_ratio+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2()); double X = Pcrit/pow(rhoc,2)*HEOS.rhomolar()/dp_drho; double tau_ref = Tc/Tref; - double dp_drho_ref = HEOS.gas_constant()*Tref*(1+2*delta*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau_ref,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau_ref,delta)); + double dp_drho_ref = HEOS.gas_constant()*Tref*(R_ratio+2*delta*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau_ref,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau_ref,delta)); double Xref = Pcrit/pow(rhoc, 2)*HEOS.rhomolar()/dp_drho_ref*Tref/HEOS.T(); num = X - Xref; @@ -756,9 +759,13 @@ long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureB double Tstar=647.096,rhostar=322,pstar=22064000,lambdastar=1e-3,mustar=1e-6; double xi; int i,j; + double R = HEOS.gas_constant()/HEOS.molar_mass(); // [J/kg/K] + double R_spec = 461.51805; //[J/kg/K] + double R_ratio = R_spec/R; Tbar = HEOS.T()/Tstar; rhobar = HEOS.keyed_output(CoolProp::iDmass)/rhostar; + double rho = HEOS.keyed_output(CoolProp::iDmass); // Dilute gas contribution lambdabar_0 = sqrt(Tbar)/(2.443221e-3+1.323095e-2/Tbar+6.770357e-3/pow(Tbar,2)-3.454586e-3/pow(Tbar,3)+4.096266e-4/pow(Tbar,4)); @@ -773,11 +780,11 @@ long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureB lambdabar_1=exp(rhobar*sum); double nu=0.630,GAMMA =177.8514,gamma=1.239,xi_0=0.13,Lambda_0=0.06,Tr_bar=1.5, - qd_bar=1/0.4,pi=3.141592654, delta = HEOS.delta(), R=461.51805;//J/kg/K + qd_bar=1/0.4,pi=3.141592654, delta = HEOS.delta(); - double drhodp = 1/(R*HEOS.T()*(1+2*rhobar*HEOS.dalphar_dDelta()+rhobar*rhobar*HEOS.d2alphar_dDelta2())); + double drhodp = 1/(R*HEOS.T()*(R_ratio+2*rhobar*HEOS.dalphar_dDelta()+rhobar*rhobar*HEOS.d2alphar_dDelta2())); double drhobar_dpbar = pstar/rhostar*drhodp; - double drhodp_Trbar = 1/(R*Tr_bar*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,1/Tr_bar,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,1/Tr_bar,delta))); + double drhodp_Trbar = 1/(R*Tr_bar*Tstar*(R_ratio+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,1/Tr_bar,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,1/Tr_bar,delta))); double drhobar_dpbar_Trbar = pstar/rhostar*drhodp_Trbar; double cp = HEOS.cpmass(); // [J/kg/K] double cv = HEOS.cvmass(); // [J/kg/K] diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 7711d52a..591f082b 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -565,7 +565,8 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE HEOS.calc_reducing_state(); const SimpleState & reduce = HEOS.get_reducing_state(); - long double R_u = HEOS.calc_gas_constant(); + long double R_u = HEOS.gas_constant(); + long double R_ratio = HEOS.calc_gas_constant_specified()/R_u; shared_ptr SatL = HEOS.SatL, SatV = HEOS.SatV; @@ -640,10 +641,10 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE long double d2alphar_ddelta2L = SatL->d2alphar_dDelta2(); long double d2alphar_ddelta2V = SatV->d2alphar_dDelta2(); - JL = deltaL * (1 + deltaL*dalphar_ddeltaL); - JV = deltaV * (1 + deltaV*dalphar_ddeltaV); - KL = deltaL*dalphar_ddeltaL + alpharL + log(deltaL); - KV = deltaV*dalphar_ddeltaV + alpharV + log(deltaV); + JL = deltaL * (R_ratio + deltaL*dalphar_ddeltaL); + JV = deltaV * (R_ratio + deltaV*dalphar_ddeltaV); + KL = deltaL*dalphar_ddeltaL + alpharL + log(deltaL)*R_ratio; + KV = deltaV*dalphar_ddeltaV + alpharV + log(deltaV)*R_ratio; PL = R_u*reduce.rhomolar*T*JL; PV = R_u*reduce.rhomolar*T*JV;