research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Replaced hard-coded values of physical constants with named constants

Browse source
This commit is contained in:
Ray Speth 2012-06-12 22:48:38 +00:00
parent aab72379fe
commit 584b9fe6f5
13 changed files with 26 additions and 101 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

18
include/cantera/thermo/DebyeHuckel.h
View file

@ -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 C<SUP>2</SUP> N<SUP>-1</SUP> m<SUP>-2</SUP>
* - e = 1.60217653E-19 C
* - F = 9.6485309E7 C kmol<SUP>-1</SUP>
* - R = 8.314472E3 kg m<SUP>2</SUP> s<SUP>-2</SUP> kmol<SUP>-1</SUP> K<SUP>-1</SUP>
* - T = 298.15 K
* - B_Debye = 3.28640E9 (kg/gmol)<SUP>1/2</SUP> m<SUP>-1</SUP>
* - \f$N_a\f$ = 6.0221415E26 kmol<SUP>-1</SUP>
*
* 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 C<SUP>2</SUP> N<SUP>-1</SUP> m<SUP>-2</SUP>
* - e = 1.60217653E-19 C
* - F = 9.6485309E7 C kmol<SUP>-1</SUP>
* - R = 8.314472E3 kg m<SUP>2</SUP> s<SUP>-2</SUP> kmol<SUP>-1</SUP> K<SUP>-1</SUP>
* - T = 298.15 K
* - B_Debye = 3.28640E9 (kg/gmol)<SUP>1/2</SUP> m<SUP>-1</SUP>
* - \f$N_a\f$ = 6.0221415E26 kmol<SUP>-1</SUP>
*
* @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;

5
include/cantera/thermo/HMWSoln.h
View file

@ -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 C<SUP>2</SUP> N<SUP>-1</SUP> m<SUP>-2</SUP>
* - e = 1.60217653E-19 C
* - F = 9.6485309E7 C kmol<SUP>-1</SUP>
* - R = 8.314472E3 kg m<SUP>2</SUP> s<SUP>-2</SUP> kmol<SUP>-1</SUP> K<SUP>-1</SUP>
* - T = 298.15 K
* - B_Debye = 3.28640E9 (kg/gmol)<SUP>1/2</SUP> m<SUP>-1</SUP>
* - \f$N_a\f$ = 6.0221415E26 kmol<SUP>-1</SUP>
*
* An example of a fixed value implementation is given below.
* @code

4
include/cantera/thermo/IdealGasPhase.h
View file

@ -103,8 +103,8 @@ namespace Cantera
* V^o_k(T,P) = \frac{R T}{P} \mbox{\quad where}
* \f]
*
* R = 8314.47215 Joules kmol<SUP>-1</SUP> K<SUP>-1</SUP>, 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 .
*
* <HR>
* <H2> Specification of Solution Thermodynamic Properties </H2>

1
include/cantera/thermo/Phase.h
View file

@ -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;

8
include/cantera/thermo/WaterProps.h
View file

@ -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)

8
samples/f77/isentropic.f
View file

@ -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

2
src/equil/vcs_VolPhase.cpp
View file

@ -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");

13
src/equil/vcs_nondim.cpp
View file

@ -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);

2
src/equil/vcs_solve.cpp
View file

@ -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),

3
src/equil/vcs_species_thermo.cpp
View file

@ -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);

6
src/equil/vcs_util.cpp
View file

@ -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",

55
src/thermo/WaterProps.cpp
View file

@ -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);

2
src/thermo/WaterPropsIAPWS.cpp
View file

@ -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) {