From 584b9fe6f550965282aabefd7458e97ac2754abc Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 12 Jun 2012 22:48:38 +0000 Subject: [PATCH] Replaced hard-coded values of physical constants with named constants --- include/cantera/thermo/DebyeHuckel.h | 18 --------- include/cantera/thermo/HMWSoln.h | 5 --- include/cantera/thermo/IdealGasPhase.h | 4 +- include/cantera/thermo/Phase.h | 1 - include/cantera/thermo/WaterProps.h | 8 +--- samples/f77/isentropic.f | 8 ++-- src/equil/vcs_VolPhase.cpp | 2 +- src/equil/vcs_nondim.cpp | 13 +++--- src/equil/vcs_solve.cpp | 2 +- src/equil/vcs_species_thermo.cpp | 3 +- src/equil/vcs_util.cpp | 6 +-- src/thermo/WaterProps.cpp | 55 +++----------------------- src/thermo/WaterPropsIAPWS.cpp | 2 +- 13 files changed, 26 insertions(+), 101 deletions(-) diff --git a/include/cantera/thermo/DebyeHuckel.h b/include/cantera/thermo/DebyeHuckel.h index 8e2edb4ff..f36f7857b 100644 --- a/include/cantera/thermo/DebyeHuckel.h +++ b/include/cantera/thermo/DebyeHuckel.h @@ -440,13 +440,8 @@ class PDSS_Water; * based on: * - \f$ \epsilon / \epsilon_0 \f$ = 78.54 * (water at 25C) - * - \f$ \epsilon_0 \f$= 8.854187817E-12 C2 N-1 m-2 - * - e = 1.60217653E-19 C - * - F = 9.6485309E7 C kmol-1 - * - R = 8.314472E3 kg m2 s-2 kmol-1 K-1 * - T = 298.15 K * - B_Debye = 3.28640E9 (kg/gmol)1/2 m-1 - * - \f$N_a\f$ = 6.0221415E26 kmol-1 * * An example of a fixed value implementation is given below. * @code @@ -1372,13 +1367,8 @@ public: * based on: * - \f$ \epsilon / \epsilon_0 \f$ = 78.54 * (water at 25C) - * - \f$ \epsilon_0 \f$= 8.854187817E-12 C2 N-1 m-2 - * - e = 1.60217653E-19 C - * - F = 9.6485309E7 C kmol-1 - * - R = 8.314472E3 kg m2 s-2 kmol-1 K-1 * - T = 298.15 K * - B_Debye = 3.28640E9 (kg/gmol)1/2 m-1 - * - \f$N_a\f$ = 6.0221415E26 kmol-1 * * @param temperature Temperature in kelvin. Defaults to -1, in which * case the temperature of the phase is assumed. @@ -1618,10 +1608,6 @@ protected: * based on: * epsilon/epsilon_0 = 78.54 * (water at 25C) - * epsilon_0 = 8.854187817E12 C2 N-1 m-2 - * e = 8.314472E3 kg m2 s-2 kmol-1 K-1 - * F = 9.6485309E7 C kmol-1 - * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 * T = 298.15 K * B_Debye = 3.28640E9 sqrt(kg/gmol)/m * @@ -1644,10 +1630,6 @@ protected: * based on: * epsilon/epsilon_0 = 78.54 * (water at 25C) - * epsilon_0 = 8.854187817E12 C2 N-1 m-2 - * e = 8.314472E3 kg m2 s-2 kmol-1 K-1 - * F = 9.6485309E7 C kmol-1 - * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 * T = 298.15 K */ double m_B_Debye; diff --git a/include/cantera/thermo/HMWSoln.h b/include/cantera/thermo/HMWSoln.h index 9cec17f31..2fca8c61a 100644 --- a/include/cantera/thermo/HMWSoln.h +++ b/include/cantera/thermo/HMWSoln.h @@ -919,13 +919,8 @@ class PDSS_Water; * based on: * - \f$ \epsilon / \epsilon_0 \f$ = 78.54 * (water at 25C) - * - \f$ \epsilon_0 \f$= 8.854187817E-12 C2 N-1 m-2 - * - e = 1.60217653E-19 C - * - F = 9.6485309E7 C kmol-1 - * - R = 8.314472E3 kg m2 s-2 kmol-1 K-1 * - T = 298.15 K * - B_Debye = 3.28640E9 (kg/gmol)1/2 m-1 - * - \f$N_a\f$ = 6.0221415E26 kmol-1 * * An example of a fixed value implementation is given below. * @code diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h index 1f38b8b37..b2c01ece0 100644 --- a/include/cantera/thermo/IdealGasPhase.h +++ b/include/cantera/thermo/IdealGasPhase.h @@ -103,8 +103,8 @@ namespace Cantera * V^o_k(T,P) = \frac{R T}{P} \mbox{\quad where} * \f] * - * R = 8314.47215 Joules kmol-1 K-1, from the 1999 CODATA convention. - * For a complete list of physical constants used within %Cantera, see \ref physConstants . + * R is the molar gas constant. For a complete list of physical constants + * used within %Cantera, see \ref physConstants . * *
*

Specification of Solution Thermodynamic Properties

diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 5208c9cb1..a0cd57f68 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -513,7 +513,6 @@ public: //! Dimensionless electrical charge of a single molecule of species k //! The charge is normalized by the the magnitude of the electron charge - //! ( \f$ e = 1.602\times 10^{-19}\f$ Coulombs). //! @param k species index doublereal charge(size_t k) const; diff --git a/include/cantera/thermo/WaterProps.h b/include/cantera/thermo/WaterProps.h index b51c2d019..1116b883c 100644 --- a/include/cantera/thermo/WaterProps.h +++ b/include/cantera/thermo/WaterProps.h @@ -223,15 +223,9 @@ public: * Nominal value at 25C and 1atm = 1.172576 sqrt(kg/gmol). * * Based on: - * epsilon/epsilon_0 = 78.54 - * (water at 25C) - * epsilon_0 = 8.854187817E-12 C2 N-1 m-2 - * e = 1.60217653E-19 C - * F = 9.6485309E7 C kmol-1 - * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 + * epsilon/epsilon_0 = 78.54 (water at 25C) * T = 298.15 K * B_Debye = 3.28640E9 sqrt(kg/gmol)/m - * Na = 6.0221415E26 * * @param T Temperature (kelvin) * @param P pressure (pascal) diff --git a/samples/f77/isentropic.f b/samples/f77/isentropic.f index 010d73369..7c2da35f1 100644 --- a/samples/f77/isentropic.f +++ b/samples/f77/isentropic.f @@ -7,7 +7,7 @@ c parameter (oneatm = 1.01325d5, NPTS = 200) double precision a(NPTS), dmach(NPTS), t(NPTS), $ ratio(NPTS) - + call newIdealGasMix('gri30.cti','gri30','') temp = 1200.d0 pres = 10.d0*oneatm @@ -52,12 +52,12 @@ c stagnation state properties double precision function soundspeed() implicit double precision (a-h,o-z) double precision meanMolarMass - parameter (R = 8314.3d0) + parameter (R = 8314.4621d0) gamma = cp_mass()/cv_mass() soundspeed = sqrt(gamma * R * temperature() $ / meanMolarMass()) return end - - + + diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 0d8e010b6..98b6de15e 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -981,7 +981,7 @@ double vcs_VolPhase::_updateVolPM() const if (m_totalMolesInert > 0.0) { if (m_gasPhase) { - double volI = m_totalMolesInert * 8314.47215 * Temp_ / Pres_; + double volI = m_totalMolesInert * Cantera::GasConstant * Temp_ / Pres_; m_totalVol += volI; } else { printf("unknown situation\n"); diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index 5b6029dac..0d1b9f5b8 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -40,13 +40,14 @@ double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const case VCS_UNITS_MKS: case VCS_UNITS_KJMOL: case VCS_UNITS_KCALMOL: - Farad = 1.602E-19 * 6.022136736e26/ (TKelvin * 8.314472E3); + Farad = Cantera::ElectronCharge * Cantera::Avogadro / + (TKelvin * Cantera::GasConstant); break; case VCS_UNITS_UNITLESS: - Farad = 1.602E-19 * 6.022136736e26; + Farad = Cantera::ElectronCharge * Cantera::Avogadro; break; case VCS_UNITS_KELVIN: - Farad = 1.602E-19 * 6.022136736e26/ (TKelvin); + Farad = Cantera::ElectronCharge * Cantera::Avogadro/ TKelvin; break; default: plogf("vcs_nondim_Farad error: unknown units: %d\n", mu_units); @@ -73,19 +74,19 @@ double VCS_SOLVE::vcs_nondimMult_TP(int mu_units, double TKelvin) const } switch (mu_units) { case VCS_UNITS_KCALMOL: - rt = TKelvin * 8.314472E-3 / 4.184; + rt = TKelvin * Cantera::GasConst_cal_mol_K * 1e-3; break; case VCS_UNITS_UNITLESS: rt = 1.0; break; case VCS_UNITS_KJMOL: - rt = TKelvin * 0.008314472; + rt = TKelvin * Cantera::GasConstant * 1e-6; break; case VCS_UNITS_KELVIN: rt = TKelvin; break; case VCS_UNITS_MKS: - rt = TKelvin * 8.314472E3; + rt = TKelvin * Cantera::GasConstant; break; default: plogf("vcs_nondimMult_TP error: unknown units: %d\n", mu_units); diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 24952e5a3..8324857e3 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -54,7 +54,7 @@ VCS_SOLVE::VCS_SOLVE() : m_totalMoleScale(1.0), m_useActCoeffJac(0), m_totalVol(0.0), - m_Faraday_dim(1.602e-19 * 6.022136736e26), + m_Faraday_dim(Cantera::ElectronCharge * Cantera::Avogadro), m_VCount(0), m_debug_print_lvl(0), m_timing_print_lvl(1), diff --git a/src/equil/vcs_species_thermo.cpp b/src/equil/vcs_species_thermo.cpp index cc4926df1..f12141e5e 100644 --- a/src/equil/vcs_species_thermo.cpp +++ b/src/equil/vcs_species_thermo.cpp @@ -231,8 +231,7 @@ VolStar_calc(size_t kglob, double TKelvin, double presPA) vol = SSStar_Vol0; break; case VCS_SSVOL_IDEALGAS: - // R J/kmol/K (2006 CODATA value) - vol= 8314.47215 * T / presPA; + vol= Cantera::GasConstant * T / presPA; break; default: plogf("%sERROR: unknown SSVol model\n", yo); diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp index db3296966..311aae7fe 100644 --- a/src/equil/vcs_util.cpp +++ b/src/equil/vcs_util.cpp @@ -679,20 +679,20 @@ double vcsUtil_gasConstant(int mu_units) double r; switch (mu_units) { case VCS_UNITS_KCALMOL: - r = 0.008314472/4.184; + r = Cantera::GasConst_cal_mol_K * 1e-3; break; case VCS_UNITS_UNITLESS: r = 1.0; break; case VCS_UNITS_KJMOL: - r = 0.008314472; + r = Cantera::GasConstant * 1e-6; break; case VCS_UNITS_KELVIN: r = 1.0; break; case VCS_UNITS_MKS: /* joules / kg-mol K = kg m2 / s2 kg-mol K */ - r = 8.314472E3; + r = Cantera::GasConstant; break; default: plogf("vcs_gasConstant error: uknown units: %d\n", diff --git a/src/thermo/WaterProps.cpp b/src/thermo/WaterProps.cpp index a0b1e7d57..d09aa3293 100644 --- a/src/thermo/WaterProps.cpp +++ b/src/thermo/WaterProps.cpp @@ -274,54 +274,9 @@ doublereal WaterProps::relEpsilon(doublereal T, doublereal P_pascal, return epsRel; } -/* - * ADebye calculates the value of A_Debye as a function - * of temperature and pressure according to relations - * that take into account the temperature and pressure - * dependence of the water density and dieletric constant. - * - * A_Debye -> this expression appears on the top of the - * ln actCoeff term in the general Debye-Huckel - * expression - * It depends on temperature. And, therefore, - * most be recalculated whenever T or P changes. - * - * A_Debye = (1/(8 Pi)) sqrt(2 Na dw / 1000) - * (e e/(epsilon R T))^3/2 - * - * Units = sqrt(kg/gmol) ~ sqrt(1/I) - * - * Nominal value = 1.172576 sqrt(kg/gmol) - * based on: - * epsilon/epsilon_0 = 78.54 - * (water at 25C) - * epsilon_0 = 8.854187817E-12 C2 N-1 m-2 - * e = 1.60217653E-19 C - * F = 9.6485309E7 C kmol-1 - * R = 8.314472E3 kg m2 s-2 kmol-1 K-1 - * T = 298.15 K - * B_Debye = 3.28640E9 sqrt(kg/gmol)/m - * Na = 6.0221415E26 - * - * ifunc = 0 return value - * ifunc = 1 return temperature derivative - * ifunc = 2 return temperature second derivative - * ifunc = 3 return pressure first derivative - * - * Verification: - * With the epsRelWater value from the BP relation, - * and the water density from the WaterDens function, - * The A_Debye computed with this function agrees with - * the Pitzer table p. 99 to 4 significant digits at 25C. - * and 20C. (Aphi = ADebye/3) - * - * (statically defined within the object) - */ + doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc) { - const doublereal e = 1.60217653E-19; - const doublereal epsilon0 = 8.854187817E-12; - const doublereal R = 8.314472E3; doublereal psat = satPressure(T); doublereal P; if (psat > P_input) { @@ -334,12 +289,12 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc) doublereal epsRelWater = relEpsilon(T, P, 0); //printf("releps calc = %g, compare to 78.38\n", epsRelWater); //doublereal B_Debye = 3.28640E9; - const doublereal Na = 6.0221415E26; - doublereal epsilon = epsilon0 * epsRelWater; + doublereal epsilon = epsilon_0 * epsRelWater; doublereal dw = density_IAPWS(T, P); - doublereal tmp = sqrt(2.0 * Na * dw / 1000.); - doublereal tmp2 = e * e * Na / (epsilon * R * T); + doublereal tmp = sqrt(2.0 * Avogadro * dw / 1000.); + doublereal tmp2 = ElectronCharge * ElectronCharge * Avogadro / + (epsilon * GasConstant * T); doublereal tmp3 = tmp2 * sqrt(tmp2); doublereal A_Debye = tmp * tmp3 / (8.0 * Pi); diff --git a/src/thermo/WaterPropsIAPWS.cpp b/src/thermo/WaterPropsIAPWS.cpp index 827d15f54..bfa3b9eb5 100644 --- a/src/thermo/WaterPropsIAPWS.cpp +++ b/src/thermo/WaterPropsIAPWS.cpp @@ -590,7 +590,7 @@ int WaterPropsIAPWS::phaseState(bool checkState) const doublereal T = T_c / tau; doublereal rho = delta * Rho_c; //doublereal psatTable = psat_est(T); - doublereal rhoMidAtm = 0.5 * (1.01E5 * M_water / (8314.472 * 373.15) + 1.0E3); + doublereal rhoMidAtm = 0.5 * (OneAtm * M_water / (Rgas * 373.15) + 1.0E3); doublereal rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15); int iStateGuess = WATER_LIQUID; if (rho < rhoMid) {