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

Fixed some errors in DebyeHuckel that weren't covered by test suite

Browse source 
(partial molar enthalpy and heat capacity)
Doxygen update
  - added private functions and added some water property files.
This commit is contained in:
Harry Moffat 2007-03-12 01:12:57 +00:00
parent 51a5f82474
commit ebdc0d014b
16 changed files with 506 additions and 188 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/src/ConstDensityThermo.h
View file

@ -443,9 +443,9 @@ namespace Cantera {
//! Current pressure (Pa)
doublereal m_press;
private:
//! Function to update the reference state thermo functions
void _updateThermo() const;
};
}

26
Cantera/src/Kinetics.cpp
View file

@ -214,19 +214,17 @@ namespace Cantera {
m_phaseindex[m_thermo.back()->id()] = nPhases();
}
/**
* err():
*
* Private function of the class Kinetics, indicating that a function
* inherited from the base class hasn't had a definition assigned to it
*/
void Kinetics::err(std::string m) const {
throw CanteraError("Kinetics::" + m,
"The default Base class method was called, when "
"the inherited class's method should "
"have been called");
}
//! Private function of the class Kinetics, indicating that a function
//! inherited from the base class hasn't had a definition assigned to it
/*!
* @param m String message
*/
void Kinetics::err(std::string m) const {
throw CanteraError("Kinetics::" + m,
"The default Base class method was called, when "
"the inherited class's method should "
"have been called");
}
}

10
Cantera/src/Kinetics.h
View file

@ -936,8 +936,14 @@ namespace Cantera {
private:
std::vector<grouplist_t> m_dummygroups;
void err(std::string m) const;
//! Vector of group lists
std::vector<grouplist_t> m_dummygroups;
//! Function for unhandled situations
/*!
* @param m String error message
*/
void err(std::string m) const;
};

14
Cantera/src/NasaThermo.h
View file

@ -468,10 +468,20 @@ namespace Cantera {
private:
//! see SpeciesThermoFactory.cpp for the definition
/*!
* @param name string name of species
* @param tmid Mid temperature, between the two temperature regions
* @param clow coefficients for lower temperature region
* @param chigh coefficients for higher temperature region
*/
void checkContinuity(std::string name, double tmid, const doublereal* clow,
doublereal* chigh);
//! for internal use by checkContinuity
/*!
* @param t temperature
* @param c coefficient array
*/
doublereal enthalpy_RT(double t, const doublereal* c) {
return c[0] + 0.5*c[1]*t + OneThird*c[2]*t*t
+ 0.25*c[3]*t*t*t + 0.2*c[4]*t*t*t*t
@ -479,6 +489,10 @@ namespace Cantera {
}
//! for internal use by checkContinuity
/*!
* @param t temperature
* @param c coefficient array
*/
doublereal entropy_R(double t, const doublereal* c) {
return c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t
+ OneThird*c[3]*t*t*t + 0.25*c[4]*t*t*t*t

91
Cantera/src/Phase.h
View file

@ -71,10 +71,10 @@ namespace Cantera {
* The XML_Node for the phase contains all of the input data used
* to set up the model for the phase, during its initialization.
*/
XML_Node& xml() { return *m_xml; }
XML_Node& xml() { return *m_xml; }
//! Return the string id for the phase
std::string id() const { return m_id; }
std::string id() const { return m_id; }
//! Set the string id for the phase
/*!
@ -92,7 +92,7 @@ namespace Cantera {
void setName(std::string nm) { m_name = nm; }
//! Returns the index of the phase
int index() const { return m_index; }
int index() const { return m_index; }
//! Sets the index of the phase
/*!
@ -329,10 +329,10 @@ namespace Cantera {
*/
doublereal massFraction(std::string name) const;
/**
* Charge density [C/m^3].
*/
doublereal chargeDensity() const;
/**
* Charge density [C/m^3].
*/
doublereal chargeDensity() const;
/// Returns the number of spatial dimensions (1, 2, or 3)
int nDim() {return m_ndim;}
@ -346,47 +346,64 @@ namespace Cantera {
*/
void setNDim(int ndim) {m_ndim = ndim;}
/**
* Finished adding species, prepare to use them for calculation
* of mixture properties.
*/
virtual void freezeSpecies();
/**
* Finished adding species, prepare to use them for calculation
* of mixture properties.
*/
virtual void freezeSpecies();
virtual bool ready() const;
virtual bool ready() const;
protected:
/**
* m_kk = Number of species in the phase. @internal m_kk is a
* member of both the State and Constituents classes.
* Therefore, to avoid multiple inheritance problems, we need
* to restate it in here, so that the declarations in the two
* base classes become hidden.
*/
int m_kk;
/**
* m_ndim is the dimensionality of the phase. Volumetric
* phases have dimensionality 3 and surface phases have
* dimensionality 2.
*/
int m_ndim;
/**
* m_index is the index of the phase
*
*/
int m_index;
/**
* m_kk = Number of species in the phase. @internal m_kk is a
* member of both the State and Constituents classes.
* Therefore, to avoid multiple inheritance problems, we need
* to restate it in here, so that the declarations in the two
* base classes become hidden.
*/
int m_kk;
/**
* m_ndim is the dimensionality of the phase. Volumetric
* phases have dimensionality 3 and surface phases have
* dimensionality 2.
*/
int m_ndim;
/**
* m_index is the index of the phase
*
*/
int m_index;
private:
vector_fp m_data;
XML_Node* m_xml;
std::string m_id;
std::string m_name;
//! This stores the initial state of the system
/*!
* @deprecated
* This doesn't seem to be used much anymore.
*/
vector_fp m_data;
//! Pointer to the XML node containing the XML info for this phase
XML_Node* m_xml;
//! ID of the phase.
/*!
* This is the value of the ID attribute of the XML phase node.
*/
std::string m_id;
//! Name of the phase.
/*!
* Initially, this is the value of the ID attribute of the XML phase node.
*/
std::string m_name;
};
//! typedef for the base Phase class
typedef Phase phase_t;
typedef Phase phase_t;
}
#endif

11
Cantera/src/ThermoPhase.h
View file

@ -1394,14 +1394,19 @@ namespace Cantera {
private:
doublereal err(std::string msg) const;
//! Error function that gets called for unhandled cases
/*!
* @param msg String containing the message.
*/
doublereal err(std::string msg) const;
};
//! typedef for the ThermoPhase class
typedef ThermoPhase thermophase_t;
typedef ThermoPhase thermophase_t;
//! typedef for the ThermoPhase class
typedef ThermoPhase thermo_t;
typedef ThermoPhase thermo_t;
}
#endif

39
Cantera/src/thermo/DebyeHuckel.cpp
View file

@ -1987,7 +1987,7 @@ namespace Cantera {
}
/**
/*
* _activityWaterHelgesonFixedForm()
*
* Formula for the log of the activity of the water
@ -2300,6 +2300,7 @@ namespace Cantera {
void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const {
double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
int k;
// First we store dAdT explicitly here
double dAdT = dA_DebyedT_TP();
if (dAdT == 0.0) {
for (k = 0; k < m_kk; k++) {
@ -2318,13 +2319,17 @@ namespace Cantera {
double sqrtI = sqrt(m_IionicMolality);
double numdAdTTmp = dAdT * sqrtI;
double denomTmp = m_B_Debye * sqrtI;
double d_lnActivitySolvent_dT = 0;
switch (m_formDH) {
case DHFORM_DILUTE_LIMIT:
for (int k = 0; k < m_kk; k++) {
for (int k = 1; k < m_kk; k++) {
m_dlnActCoeffMolaldT[k] =
m_lnActCoeffMolal[k] * dAdT / m_A_Debye;
}
d_lnActivitySolvent_dT = 2.0 / 3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrt(m_IionicMolality);
m_dlnActCoeffMolaldT[m_indexSolvent] = d_lnActivitySolvent_dT;
break;
case DHFORM_BDOT_AK:
@ -2365,7 +2370,8 @@ namespace Cantera {
sigma = 0.0;
}
m_dlnActCoeffMolaldT[m_indexSolvent] =
2.0 /3.0 * dAdT * m_Mnaught * m_IionicMolality * sqrtI * sigma;
2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
@ -2384,8 +2390,7 @@ namespace Cantera {
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldT[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
m_dlnActCoeffMolaldT[m_indexSolvent] =
2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@ -2405,8 +2410,7 @@ namespace Cantera {
}
sigma = 1.0 / ( 1.0 + denomTmp);
m_dlnActCoeffMolaldT[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
m_dlnActCoeffMolaldT[m_indexSolvent] =
2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@ -2416,6 +2420,8 @@ namespace Cantera {
exit(-1);
break;
}
}
/*
@ -2521,8 +2527,7 @@ namespace Cantera {
} else {
sigma = 0.0;
}
m_d2lnActCoeffMolaldT2[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@ -2533,17 +2538,16 @@ namespace Cantera {
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldT[k] =
m_d2lnActCoeffMolaldT2[k] =
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * d2AdT2 * tmpLn
/ (m_B_Debye * m_Aionic[0]);
m_dlnActCoeffMolaldT[k] /= 3.0;
m_d2lnActCoeffMolaldT2[k] /= 3.0;
}
}
sigma = 1.0 / ( 1.0 + denomTmp);
m_dlnActCoeffMolaldT[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@ -2641,7 +2645,8 @@ namespace Cantera {
sigma = 0.0;
}
m_dlnActCoeffMolaldP[m_indexSolvent] =
2.0 /3.0 * dAdP * m_Mnaught * m_IionicMolality * sqrtI * sigma;
2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
@ -2660,8 +2665,7 @@ namespace Cantera {
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldP[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
m_dlnActCoeffMolaldP[m_indexSolvent] =
2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@ -2681,8 +2685,7 @@ namespace Cantera {
}
sigma = 1.0 / ( 1.0 + denomTmp);
m_dlnActCoeffMolaldP[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
m_dlnActCoeffMolaldP[m_indexSolvent] =
2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;

41
Cantera/src/thermo/DebyeHuckel.h
View file

@ -53,7 +53,9 @@ namespace Cantera {
* The concentrations of the ionic species are assumed to obey the electroneutrality
* condition.
*
* <b> Specification of Species Standard %State Properties </b>
* <HR>
* <H2> Specification of Species Standard %State Properties </H2>
* <HR>
*
* The standard states are on the unit molality basis. Therefore, in the
* documentation below, the normal \f$ o \f$ superscript is replaced with
@ -454,10 +456,10 @@ namespace Cantera {
* Currently, \f$ B_{Debye} \f$ is a constant in the model, specified either by a default
* water value, or through the input file. This may have to be looked at, in the future.
*
*
* <HR>
* <H2> %Application within %Kinetics Managers </H2>
* <HR>
*
* For the time being, we have set the standard concentration for all species in
* this phase equal to the default concentration of the solvent at 298 K and 1 atm.
* This means that the
@ -507,7 +509,7 @@ namespace Cantera {
* Note, this treatment may be modified in the future, as events dictate.
*
* <HR>
* <b> Instantiation of the Class </b>
* <H2> Instantiation of the Class </H2>
* <HR>
*
* The constructor for this phase is NOT located in the default ThermoFactory
@ -540,7 +542,7 @@ namespace Cantera {
* @endcode
*
* <HR>
* <b> XML Example </b>
* <H2> XML Example </H2>
* <HR>
*
* The phase model name for this is called StoichSubstance. It must be supplied
@ -597,8 +599,6 @@ namespace Cantera {
</phase>
@endverbatim
*
* The model attribute, "StoichSubstanceSSTP", on the thermo element identifies the phase as
* being a StoichSubstanceSSTP object.
*
*/
class DebyeHuckel : public MolalityVPSSTP {
@ -994,16 +994,19 @@ namespace Cantera {
//! in the mixture. Units (J/kmol)
/*!
* For this phase, the partial molar enthalpies are equal to the
* pure species enthalpies
* standard state enthalpies modified by the derivative of the
* molality-based activity coefficent wrt temperature
*
* \f[
* \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k
* \bar h_k(T,P) = h^{\triangle}_k(T,P) - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
* \f]
* The reference-state pure-species enthalpies,
* \f$ \hat h^{ref}_k(T) \f$,
* at the reference pressure,\f$ P_{ref} \f$,
* are computed by the species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
* The solvent partial molar enthalpy is equal to
* \f[
* \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o}{dT}
* \f]
*
* The temperature dependence of the activity coefficients currently
* only occurs through the temperature dependence of the Debye constant.
*
* @param hbar Output vector of species partial molar enthalpies.
* Length: m_kk. units are J/kmol.
@ -1562,6 +1565,8 @@ namespace Cantera {
//! salt-out modifications.
/*!
* Returns the calculated activity coefficients.
*
* @param IionicMolality Value of the ionic molality (sqrt(gmol/kg))
*/
double _nonpolarActCoeff(double IionicMolality) const;
@ -1572,6 +1577,12 @@ namespace Cantera {
* NaCl brine. It's to be used with extreme caution.
*/
double _osmoticCoeffHelgesonFixedForm() const;
//! Formula for the log of the water activity that occurs in the GWB.
/*!
* It is originally from Helgeson for a variable
* NaCl brine. It's to be used with extreme caution.
*/
double _lnactivityWaterHelgesonFixedForm() const;
@ -1868,6 +1879,8 @@ namespace Cantera {
*
* We assume that the activity coefficients are current in this routine
*
*
*
* The solvent activity coefficient is on the molality scale. It's derivative is too.
*/
void s_update_dlnMolalityActCoeff_dT() const;

6
Cantera/src/thermo/SingleSpeciesTP.h
View file

@ -703,8 +703,12 @@ namespace Cantera {
void _updateThermo() const;
private:
doublereal err(std::string msg) const;
//! Error return for unhandled cases
/*!
* @param msg String message
*/
doublereal err(std::string msg) const;
};
}

6
Cantera/src/thermo/VPStandardStateTP.h
View file

@ -572,12 +572,12 @@ namespace Cantera {
*/
mutable vector_fp m_Vss;
private:
//! VPStandardStateTP has its own err routine
/*!
* VPStandardStateTP has its own err routine
* @param msg Error message string
*/
doublereal err(std::string msg) const;

10
Cantera/src/thermo/WaterPropsIAPWS.cpp
View file

@ -78,7 +78,7 @@ double WaterPropsIAPWS::helmholtzFE_RT() const{
return (retn);
}
/**
/*
* Calculate the Helmholtz free energy in mks units of
* J kmol-1 K-1.
*/
@ -96,7 +96,7 @@ double WaterPropsIAPWS::helmholtzFE() const{
}
/**
/*
* Calculate the pressure (Pascals), given the temperature and density
* Temperature: kelvin
* rho: density in kg m-3
@ -113,7 +113,7 @@ double WaterPropsIAPWS::pressure() const{
return (retn * rho * Rgas * temperature);
}
/**
/*
* Calculates the pressure in dimensionless form
* p/(rhoRT) at the currently stored tau and delta values
*/
@ -199,7 +199,7 @@ double WaterPropsIAPWS::density() const {
return (delta * Rho_c);
}
/**
/*
* psat_est provides a rough estimate of the saturation
* pressure given the temperature. This is used as an initial
* guess for refining the pressure.
@ -245,7 +245,7 @@ double WaterPropsIAPWS::psat_est(double temperature) {
return ps;
}
/**
/*
* Returns the coefficient of thermal expansion as a function
* of temperature and pressure.
* alpha = d (ln V) / dT at constant P.

284
Cantera/src/thermo/WaterPropsIAPWS.h
View file

@ -16,15 +16,19 @@
#include "WaterPropsIAPWSphi.h"
/*
* These constants are defined and used in the interphase
* to describe desired phases.
*/
/**
* @name Names for the phase regions
*
* These constants are defined and used in the interface
* to describe desired phases.
*/
//@{
#define WATER_GAS 0
#define WATER_LIQUID 1
#define WATER_SUPERCRIT 2
//@}
/**
/*!
* Class for calculating the properties of water.
*
*
@ -32,149 +36,264 @@
* used in the steam tables, i.e., the liquid at the triple point
* for water has the following properties:
*
* u(273.16, rho) = 0.0
* s(273.16, rho) = 0.0
* psat(273.16) = 611.655 Pascal
* rho(273.16, psat) = 999.793 kg m-3
* - u(273.16, rho) = 0.0
* - s(273.16, rho) = 0.0
* - psat(273.16) = 611.655 Pascal
* - rho(273.16, psat) = 999.793 kg m-3
*
*/
class WaterPropsIAPWS {
public:
//! Base constructor
WaterPropsIAPWS();
WaterPropsIAPWS(const WaterPropsIAPWS &b);
WaterPropsIAPWS & operator=(const WaterPropsIAPWS &b);
//! Copy constructor
WaterPropsIAPWS(const WaterPropsIAPWS &);
//! assignment constructor
WaterPropsIAPWS & operator=(const WaterPropsIAPWS &);
//! destructor
~WaterPropsIAPWS();
//! Set the internal state of the object wrt temperature and density
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
void setState(double temperature, double rho);
/**
* Calculate the Helmholtz free energy in mks units of
* J kmol-1 K-1.
//! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1.
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double helmholtzFE(double temperature, double rho);
//! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1,
//! using the last temperature and density
double helmholtzFE() const;
/**
* Calculate the Gibbs free energy in mks units of
* J kmol-1 K-1.
//! Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double Gibbs(double temperature, double rho);
//! Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
//! using the last temperature and density
double Gibbs() const;
/**
* Calculate the enthalpy in mks units of
* J kmol-1
//! Calculate the enthalpy in mks units of J kmol-1
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double enthalpy(double temperature, double rho);
//! Calculate the enthalpy in mks units of J kmol-1
//! using the last temperature and density
double enthalpy() const;
/**
* Calculate the internal energy in mks units of
* J kmol-1
//! Calculate the internal energy in mks units of J kmol-1
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double intEnergy(double temperature, double rho);
double intEnergy() const;
/**
* Calculate the entropy in mks units of
* J kmol-1 K-1
//! Calculate the internal energy in mks units of J kmol-1
//! at the last internal energy
double intEnergy() const;
//! Calculate the entropy in mks units of J kmol-1 K-1
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double entropy(double temperature, double rho);
//! Calculate the entropy in mks units of J kmol-1 K-1
//! at the last temperature and density
double entropy() const;
/**
* Calculate the constant volume heat capacity
* in mks units of J kmol-1 K-1
//! Calculate the constant volume heat capacity in mks units of J kmol-1 K-1
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double cv(double temperature, double rho);
//! Calculate the constant volume heat capacity in mks units of J kmol-1 K-1
//! at the last temperature and density
double cv() const;
/**
* Calculate the constant pressure heat capacity
* in mks units of J kmol-1 K-1
//! Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double cp(double temperature, double rho);
//! Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1
//! at the last temperature and density
double cp() const;
//! Calculate the molar volume (kmol m-3)
/*!
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
double molarVolume(double temperature, double rho);
double molarVolume() const;
/**
* Calculate the pressure (Pascals), given the temperature and density
* Temperature: kelvin
* rho: density in kg m-3
//! Calculate the molar volume (kmol m-3)
//! at the last temperature and density
double molarVolume() const;
//! Calculate the pressure (Pascals), given the temperature and density
/*!
* @param temperature input temperature kelvin
* @param rho density in kg m-3
*
* @return
* returns the pressure (Pascal)
*/
double pressure(double temperature, double rho);
//! Calculates the pressure (Pascals), given the current value of the
//! temperature and density.
/*!
* The density is an independent variable in the underlying equation of state
*/
double pressure() const;
/*
* Calculates the density given the temperature and the pressure,
* and a guess at the density. Note, below T_c, this is a
* multivalued function.
//! Calculates the density given the temperature and the pressure,
//! and a guess at the density.
/*!
* Note, below T_c, this is a multivalued function.
*
* parameters:
* temperature: Kelvin
* pressure : Pressure in Pascals (Newton/m**2)
* phase : guessed phase of water
* : -1: no guessed phase
* rhoguess : guessed density of the water
* : -1.0 no guessed density
*
* @param temperature: Kelvin
* @param pressure : Pressure in Pascals (Newton/m**2)
* @param phase : guessed phase of water
* : -1: no guessed phase
* @param rhoguess : guessed density of the water
* : -1.0 no guessed density
* @return
* Returns the density
*/
double density(double temperature, double pressure,
int phase = -1, double rhoguess = -1.0);
//! Returns the density (kg m-3)
/*!
* The density is an independent variable in the underlying equation of state
*/
double density() const;
/**
* This function returns an estimated value for the saturation
* pressure. It does this via a polynomial fit of the vapor pressure
* curve.
//! This function returns an estimated value for the saturation pressure.
/*!
* It does this via a polynomial fit of the vapor pressure curve.
* units = (Pascals)
*
* @param temperature Input temperature (Kelvin)
*
* @return
* Returns the estimated saturation pressure
*/
double psat_est(double temperature);
/**
* Returns the coefficient of thermal expansion as a function
* of temperature and pressure.
//! Returns the coefficient of thermal expansion as a function of temperature and pressure.
/*!
* alpha = d (ln V) / dT at constant P.
*
*
* @param temperature Input temperature (Kelvin)
* @param pressure Input pressure (Pa)
* @return
* Returns the coefficient of thermal expansion
*/
double coeffThermExp(double temperature, double pressure);
/**
* Returns the coefficient of isothermal compressibility as a function
* of temperature and pressure.
//! Returns the coefficient of isothermal compressibility as a function
//! of temperature and pressure.
/*!
* kappa = - d (ln V) / dP at constant T.
*
* units - 1/Pascal
*
* @param temperature Input temperature (Kelvin)
* @param pressure Input pressure (Pa)
* @return
* returns the isothermal compressibility
*/
double isothermalCompressibility(double temperature, double pressure);
/**
* Utility routine in the calculation of the saturation pressure
*/
//! Utility routine in the calculation of the saturation pressure
/*!
* @param temperature temperature (kelvin)
* @param pressure pressure (Pascal)
* @param densLiq Output density of liquid
* @param densGas output Density of gas
* @param delGRT output delGRT
*/
void corr(double temperature, double pressure, double &densLiq,
double &densGas, double &delGRT);
//! Utility routine in the calculation of the saturation pressure
/*!
* @param temperature temperature (kelvin)
* @param pressure pressure (Pascal)
* @param densLiq Output density of liquid
* @param densGas output Density of gas
* @param pcorr output corrected pressure
*/
void corr1(double temperature, double pressure, double &densLiq,
double &densGas, double &pcorr);
/**
* This function returns the saturation pressure given the
* temperature as an input parameter.
* units = Pascal
//! This function returns the saturation pressure given the
//! temperature as an input parameter.
/*!
* @param temperature input temperature (kelvin)
* @return
* Returns the saturation pressure
* units = Pascal
*/
double psat(double temperature);
//! Returns the critical temperature of water (Kelvin)
/*!
* This is hard coded to the value 647.096 Kelvin
*/
double Tcrit() { return 647.096;}
//! Returns the critical pressure of water (22.064E6 Pa)
/*!
* This is hard coded to the value of 22.064E6 pascals
*/
double Pcrit() { return 22.064E6;}
//! Return the critical density of water (kg m-3)
/*!
* This is equal to 322 kg m-3.
*/
double Rhocrit() { return 322.;}
private:
/**
* Calculate the dimensionless temp and rho and store internally.
*
* @param temperature input temperature (kelvin)
* @param rho density in kg m-3
*/
void calcDim(double temperature, double rho);
@ -184,19 +303,46 @@ private:
* show the dimensional functions in the interface.
*/
double helmholtzFE_RT() const;
//! Returns the dimensionless gibbs free energy
double Gibbs_RT() const;
//! Returns the dimensionless enthalpy
double enthalpy_RT() const;
//! Returns the dimensionless internal energy
double intEnergy_RT() const;
//! Returns the dimensionless entropy
double entropy_R() const;
//! Returns the dimensionless heat capacity at constant volume
double cv_R() const;
//! Returns the dimensionless heat capacity at constant pressure
double cp_R() const;
//! Return the current dimensionless pressure
double pressure_rhoRT() const;
protected:
//! pointer to the underlying object that does the calculations.
WaterPropsIAPWSphi *m_phi;
//! Dimensionless temperature
/*!
* tau = T_C / T
*/
double tau;
//! Dimensionless density
/*!
* delta = rho / rho_c
*/
double delta;
//! Current state of the system
int iState;
};
#endif

18
Cantera/src/thermo/WaterPropsIAPWSphi.cpp
View file

@ -363,7 +363,7 @@ WaterPropsIAPWSphi::WaterPropsIAPWSphi() :
{
}
/**
/*
* intCheck() calculates all of the functions at a one point and
* prints out the result. It's used for conducting the internal
* check.
@ -624,7 +624,7 @@ double WaterPropsIAPWSphi::phi_d(double tau, double delta) {
return retn;
}
/**
/*
* Calculate the dimensionless pressure at tau and delta;
*
* p/(rhoRT) = delta * phi_d()
@ -854,7 +854,7 @@ double WaterPropsIAPWSphi::phiR_t() const {
return val;
}
/**
/*
* Calculate the dPhidtau function, which is basically the derivative
* of helmholtz free energy wrt tau
* Eqn. (6.4)
@ -867,7 +867,7 @@ double WaterPropsIAPWSphi::phi_t(double tau, double delta) {
return retn;
}
/**
/*
* Calculate d2_phi0/dtau2
*/
double WaterPropsIAPWSphi::phi0_tt() const {
@ -882,7 +882,7 @@ double WaterPropsIAPWSphi::phi0_tt() const {
return retn;
}
/**
/*
* Calculate Eqn. 6.6 for dphiRdtau, the second derivative residual part of the
* dimensionless Helmholtz free energy wrt temperature
*
@ -1190,7 +1190,7 @@ double WaterPropsIAPWSphi::enthalpy_RT() const {
return hRT;
}
/**
/*
* Calculate the dimensionless entropy s/R.
*/
double WaterPropsIAPWSphi::entropy_R() const {
@ -1203,7 +1203,7 @@ double WaterPropsIAPWSphi::entropy_R() const {
return sR;
}
/**
/*
* Calculate the dimensionless internal energy, u/RT.
*/
double WaterPropsIAPWSphi::intEnergy_RT() const {
@ -1214,7 +1214,7 @@ double WaterPropsIAPWSphi::intEnergy_RT() const {
return uR;
}
/**
/*
* Calculate the dimensionless constant volume Heat Capacity, Cv/R
*/
double WaterPropsIAPWSphi::cv_R() const {
@ -1225,7 +1225,7 @@ double WaterPropsIAPWSphi::cv_R() const {
return cvR;
}
/**
/*
* Calculate the dimensionless constant pressure Heat Capacity, Cp/R
*/
double WaterPropsIAPWSphi::cp_R() const {

111
Cantera/src/thermo/WaterPropsIAPWSphi.h
View file

@ -1,5 +1,7 @@
/**
* @file WaterPropsIAPWSphi.h
*
* Lowest level of the classes which support a real water model.
*/
/*
* Copywrite (2006) Sandia Corporation. Under the terms of
@ -13,31 +15,81 @@
#ifndef WATERPROPSIAPWSPHI_H
#define WATERPROPSIAPWSPHI_H
/*
/*!
* the WaterPropsIAPSWSphi class support low level calls for
* the real description of water.
*
* Units Note: This class works with reduced units exclusively.
*/
class WaterPropsIAPWSphi {
public:
//! Base constructor
WaterPropsIAPWSphi();
/*
* Calculate the base phi's, recalculating the internal polynomials
//! Calculate the Phi function, which is the base function
/*!
* The phi functino is basically the helmholtz free energy
* Eqn. (6.4)
* All internal polynomials are recalculated.
*
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
double phi(double tau, double delta);
double phi_d(double tau, double delta);
//! Delta derivative of phi
/*!
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
double phi_d(double tau, double delta);
//! 2nd derivative of phi wrt delta
/*!
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
double phi_dd(double tau, double delta);
//! First derivative of phi wrt tau
/*!
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
double phi_t(double tau, double delta);
//! Second derivative of phi wrt tau
/*!
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
double phi_tt(double tau, double delta);
//! Second derivative of phi wrt tau, then delta
/*!
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
double phi_dt(double tau, double delta);
//! Internal check # 1
void check1();
//! Internal check # 2
void check2();
/**
* Calculate the dimensionless pressure, pred:
* pred = pressure M / (rho RT)
//! Calculate the dimensionless pressure at tau and delta;
/*!
*
* p/(rhoRT) = delta * phi_d() = 1.0 + delta phiR_d()
*
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*
* note: this is done so much, we have a seperate routine.
*/
double pressure_rhoRT(double tau, double delta);
@ -46,6 +98,13 @@ public:
* and the reduced temperature, tau. It takes an initial guess, deltaGuess.
* DeltaGuess is important as this is a multivalued function below the
* critical point.
*
* @param p_red Value of the dimensionless pressure
* @param tau Dimensionless temperature = T_c/T
* @param deltaGuess Initial guess for the dimensionless density
*
* @return
* Returns the dimensionless density.
*/
double dfind(double p_red, double tau, double deltaGuess);
@ -83,29 +142,65 @@ public:
* Calculates internal polynomials in tau and delta. This
* routine is used to store the internal state of tau and delta
* for later use by the other routines in the class.
*
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
void tdpolycalc(double tau, double delta);
//! Return the value of phiR(), res
double phiR() const;
private:
//! nau calculation
double phi0() const;
//! calculation of d_phiR/d_d
double phiR_d() const;
//! calculation of d_nau/d_d
double phi0_d() const;
//! calculation of d2_res/d_dd
double phiR_dd() const;
//! calculation of d2_nau/d_dd
double phi0_dd() const;
//! calculation of d_nau/d_t
double phi0_t() const;
//! calculation of d_res/d_t
double phiR_t() const;
//! calculation of d2_res/d_tt
double phiR_tt() const;
//! calculation of d2_nau/d_tt
double phi0_tt() const;
//! calculation of d2_res/d_dt
double phiR_dt() const;
//! calculation of d2_nau/d_dt
double phi0_dt() const;
/**
* intCheck() calculates all of the functions at a one point and
* prints out the result. It's used for conducting the internal
* check.
*
* @param tau Dimensionless temperature = T_c/T
* @param delta Dimensionless density = delta = rho / Rho_c
*/
void intCheck(double tau, double delta);
protected:
//! Value of internally calculated polynomials of powers of TAU
double TAUp[52];
//! Value of internally calculated polynomials of powers of delta
double DELTAp[16];
//! Last tau that was used to calculate polynomials
double TAUsave;
//! sqrt of TAU
double TAUsqrt;
//! Last delta that was used to calculate polynomials
double DELTAsave;
};
#endif

20
Cantera/src/units.h
View file

@ -129,9 +129,25 @@ namespace Cantera {
/// pointer to the single instance of Unit
static Unit* s_u;
std::map<std::string, doublereal> m_u;
std::map<std::string, doublereal> m_act_u;
//! Map between a string and a units double value
/*!
* This map maps the dimension string to the units value adjustment. Example
* - m_u["m"] = 1.0;
* - m_u["cm"] = 0.01;
*/
std::map<std::string, doublereal> m_u;
//! Map between a string and a units double value for activation energy units
/*!
* This map maps the dimension string to the units value adjustment. Example
* - m_act_u["K"] = GasConstant;
*/
std::map<std::string, doublereal> m_act_u;
/*!
* Units class constructor, containing the default mappings between
* strings and units.
*/
Unit(){
// length

5
tools/doc/Cantera.cfg.in
View file

@ -43,7 +43,7 @@ SUBGROUPING = YES
# Build related configuration options
#---------------------------------------------------------------------------
EXTRACT_ALL = NO
EXTRACT_PRIVATE = NO
EXTRACT_PRIVATE = YES
EXTRACT_STATIC = NO
EXTRACT_LOCAL_CLASSES = YES
EXTRACT_LOCAL_METHODS = NO
@ -117,7 +117,8 @@ FILE_PATTERNS = Kinetics.h Kinetics.cpp \
IdealSolidSolnPhase.h IdealSolidSolnPhase.cpp \
StoichSubstanceSSTP.h StoichSubstanceSSTP.cpp \
DebyeHuckel.h DebyeHuckel.cpp \
ConstDensityThermo.h ConstDensityThermo.cpp
ConstDensityThermo.h ConstDensityThermo.cpp \
WaterPropsIAPWSphi.h WaterPropsIAPWSphi.cpp WaterPropsIAPWS.h WaterPropsIAPWS.cpp
RECURSIVE = NO
EXCLUDE = CVS examples converters zeroD
EXCLUDE_SYMLINKS = NO