Incremental update on HKFT implementation

Added thermalExpansionCoeff() and its temp derivative to objects.
This commit is contained in:
Harry Moffat 2008-08-21 01:46:11 +00:00
parent 045aa04116
commit 96930d57d6
8 changed files with 293 additions and 158 deletions

View file

@ -356,11 +356,13 @@ namespace Cantera {
double t2 = temp * temp;
double val = ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
}
if (ifunc != 1) {
throw CanteraError("HKFT_PDSS::ag", "unimplemented");
if (ifunc != 2) {
return 0.0;
}
return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
return ag_coeff[2] * 2.0;;
}
//! Internal formula for the calculation of b_g()
@ -373,11 +375,13 @@ namespace Cantera {
double t2 = temp * temp;
double val = bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2;
return val;
} else if (ifunc == 1) {
return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
}
if (ifunc != 1) {
throw CanteraError("HKFT_PDSS::bg", "unimplemented");
if (ifunc != 2) {
return 0.0;
}
return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
return bg_coeff[2] * 2.0;
}
double HKFT_PDSS::f(const double temp, const double pres, const int ifunc) {
@ -402,7 +406,14 @@ namespace Cantera {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
return fac1 * fac2;
} else if (ifunc == 1) {
fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 16.0)) / 300.;
fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300.;
return fac1 * fac2;
} else if (ifunc == 2) {
fac1 = (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.);
return fac1 * fac2;
} else if (ifunc == 3) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
fac2 = (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3 )/ 1.0E5;
return fac1 * fac2;
} else {
throw CanteraError("HKFT_PDSS::gg", "unimplemented");
@ -417,18 +428,58 @@ namespace Cantera {
m_densWaterSS = m_waterSS->density();
// density in gm cm-3
double dens = m_densWaterSS * 1.0E-3;
if (ifunc == 0) {
if (dens >= 1.0) {
return 0.0;
double gval = afunc * pow((1.0-dens), bfunc);
if (dens >= 1.0) {
return 0.0;
}
if (ifunc == 0) {
return gval;
} else if (ifunc == 1 || ifunc == 2) {
double afuncdT = ag(temp, 1);
double bfuncdT = bg(temp, 1);
double alpha = m_waterSS->thermalExpansionCoeff();
double fac1 = afuncdT * gval / afunc;
double fac2 = bfuncdT * gval * log(1.0 - dens);
double fac3 = gval * alpha * bfunc * dens / (1.0 - dens);
double dgdt = fac1 + fac2 + fac3;
if (ifunc == 1) {
return dgdt;
}
double afuncdT2 = ag(temp, 2);
double bfuncdT2 = bg(temp, 2);
double dfac1dT = dgdt * afuncdT / afunc + afuncdT2 * gval / afunc
- afuncdT * afuncdT * gval / (afunc * afunc);
double ddensdT = - alpha * dens;
double dfac2dT = bfuncdT2 * gval * log(1.0 - dens)
+ bfuncdT * dgdt * log(1.0 - dens)
- bfuncdT * gval /(1.0 - dens) * ddensdT;
double dalphadT = m_waterSS->dthermalExpansionCoeffdT();
double gval = afunc * pow((1.0-dens), bfunc);
double fval = f(temp, pres, ifunc);
return gval - fval;
double dfac3dT = dgdt * alpha * bfunc * dens / (1.0 - dens)
+ gval * dalphadT * bfunc * dens / (1.0 - dens)
+ gval * alpha * bfuncdT * dens / (1.0 - dens)
+ gval * alpha * bfunc * ddensdT / (1.0 - dens)
- gval * alpha * bfunc * dens / ((1.0 - dens) * (1.0 - dens)) * ddensdT;
return dfac1dT + dfac2dT + dfac3dT;
} else {
throw CanteraError("HKFT_PDSS::gg", "unimplemented");
}
return 0.0;
}
double HKFT_PDSS::gstar(const double temp, const double pres, const int ifunc) {
double gval = g(temp, pres, ifunc);
double fval = f(temp, pres, ifunc);
return gval - fval;
}
}

View file

@ -144,6 +144,7 @@ namespace Cantera {
double bg(const double temp, const int ifunc = 0) const;
double g(const double temp, const double pres, const int ifunc = 0);
double f(const double temp, const double pres, const int ifunc = 0);
double gstar(const double temp, const double pres, const int ifunc = 0);
protected:

View file

@ -397,7 +397,18 @@ namespace Cantera {
throw CanteraError("PDSS::pressure()", "unimplemented");
}
// Return the volumetric thermal expansion coefficient. Units: 1/K.
/*
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
doublereal PDSS::thermalExpansionCoeff() const {
throw CanteraError("PDSS::thermalExpansionCoeff()", "unimplemented");
return (0.0);
}
/// critical temperature
doublereal PDSS::critTemperature() const {
throw CanteraError("PDSS::critTemperature()", "unimplemented");

View file

@ -213,6 +213,15 @@ namespace Cantera {
virtual doublereal pressure() const;
virtual void setPressure(doublereal p);
//! Return the volumetric thermal expansion coefficient. Units: 1/K.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal thermalExpansionCoeff() const;
//@}
/// @name Partial Molar Properties of the Solution -----------------
//@{

View file

@ -387,6 +387,28 @@ namespace Cantera {
}
// Return the volumetric thermal expansion coefficient. Units: 1/K.
/*
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
doublereal WaterPDSS::thermalExpansionCoeff() const {
doublereal pres = pressure();
doublereal val = m_sub->coeffThermExp(m_temp, pres);
return val;
}
doublereal WaterPDSS::dthermalExpansionCoeffdT() const {
doublereal pres = pressure();
double tt = m_temp - 0.04;
doublereal vald = m_sub->coeffThermExp(tt, pres);
doublereal val2 = m_sub->coeffThermExp(m_temp, pres);
doublereal val = (val2 - vald) / 0.04;
return val;
}
/// critical temperature
doublereal WaterPDSS::critTemperature() const { return m_sub->Tcrit(); }

View file

@ -25,178 +25,196 @@ class WaterPropsIAPWS;
namespace Cantera {
/**
* Class for the liquid water pressure dependent
* standard state
*
*
* Notes:
* Base state for thermodynamic properties:
*
* The thermodynamic base state for water is set to the NIST basis here
* by specifying constants EW_Offset and SW_Offset. These offsets are
* specified so that the following properties hold:
*
* Delta_Hfo_gas(298.15) = -241.826 kJ/gmol
* So_gas(298.15, 1bar) = 188.835 J/gmolK
*
* (http://webbook.nist.gov)
*
* The "o" here refers to a hypothetical ideal gas state. The way
* we achieve this in practice is to evaluate at a very low pressure
* and then use the theoretical ideal gas results to scale up to
* higher pressures:
*
* Ho(1bar) = H(P0)
*
* So(1bar) = S(P0) + RT ln(1bar/P0)
*
* The offsets used in the steam tables are different than NIST's.
* They assume u_liq(TP) = 0.0, s_liq(TP) = 0.0, where TP is the
* triple point conditions.
*
*
*/
class WaterPDSS : public PDSS {
public:
/**
* Class for the liquid water pressure dependent
* standard state
*
*
* Notes:
* Base state for thermodynamic properties:
*
* The thermodynamic base state for water is set to the NIST basis here
* by specifying constants EW_Offset and SW_Offset. These offsets are
* specified so that the following properties hold:
*
* Delta_Hfo_gas(298.15) = -241.826 kJ/gmol
* So_gas(298.15, 1bar) = 188.835 J/gmolK
*
* (http://webbook.nist.gov)
*
* The "o" here refers to a hypothetical ideal gas state. The way
* we achieve this in practice is to evaluate at a very low pressure
* and then use the theoretical ideal gas results to scale up to
* higher pressures:
*
* Ho(1bar) = H(P0)
*
* So(1bar) = S(P0) + RT ln(1bar/P0)
*
* The offsets used in the steam tables are different than NIST's.
* They assume u_liq(TP) = 0.0, s_liq(TP) = 0.0, where TP is the
* triple point conditions.
*
*
* Basic list of constructors and duplicators
*/
class WaterPDSS : public PDSS {
public:
/**
* Basic list of constructors and duplicators
*/
WaterPDSS();
WaterPDSS(ThermoPhase *tp, int spindex);
WaterPDSS(const WaterPDSS &b);
WaterPDSS& operator=(const WaterPDSS&b);
WaterPDSS(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id = "");
WaterPDSS(ThermoPhase *tp, int spindex,
XML_Node& phaseRef, std::string id = "");
virtual ~WaterPDSS();
WaterPDSS();
WaterPDSS(ThermoPhase *tp, int spindex);
WaterPDSS(const WaterPDSS &b);
WaterPDSS& operator=(const WaterPDSS&b);
WaterPDSS(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id = "");
WaterPDSS(ThermoPhase *tp, int spindex,
XML_Node& phaseRef, std::string id = "");
virtual ~WaterPDSS();
/**
*
* @name Utilities
* @{
*/
virtual int pdssType() const { return -1; }
/**
*
* @name Utilities
* @{
*/
virtual int pdssType() const { return -1; }
/**
* @}
* @name Molar Thermodynamic Properties of the Solution --------------
* @{
*/
virtual doublereal enthalpy_mole() const;
virtual doublereal intEnergy_mole() const;
virtual doublereal entropy_mole() const;
virtual doublereal gibbs_mole() const;
virtual doublereal cp_mole() const;
virtual doublereal cv_mole() const;
/**
* @}
* @name Molar Thermodynamic Properties of the Solution --------------
* @{
*/
virtual doublereal enthalpy_mole() const;
virtual doublereal intEnergy_mole() const;
virtual doublereal entropy_mole() const;
virtual doublereal gibbs_mole() const;
virtual doublereal cp_mole() const;
virtual doublereal cv_mole() const;
//@}
/// @name Mechanical Equation of State Properties ---------------------
//@{
//@}
/// @name Mechanical Equation of State Properties ---------------------
//@{
virtual doublereal pressure() const;
virtual void setTempPressure(doublereal t, doublereal p);
virtual void setPressure(doublereal p);
virtual doublereal pressure() const;
virtual void setTempPressure(doublereal t, doublereal p);
virtual void setPressure(doublereal p);
//@}
/// @name Partial Molar Properties of the Solution -----------------
//@{
//! Return the volumetric thermal expansion coefficient. Units: 1/K.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal thermalExpansionCoeff() const;
//! Return the derivative of the volumetric thermal expansion coefficient. Units: 1/K2.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal dthermalExpansionCoeffdT() const;
virtual void getChemPotentials(doublereal* mu) const {
mu[0] = gibbs_mole();
}
//@}
/// @name Partial Molar Properties of the Solution -----------------
//@{
//@}
/// @name Properties of the Standard State of the Species
// in the Solution --
//@{
virtual void getChemPotentials(doublereal* mu) const {
mu[0] = gibbs_mole();
}
//@}
/// @name Properties of the Standard State of the Species
// in the Solution --
//@{
/// critical temperature
virtual doublereal critTemperature() const;
/// critical temperature
virtual doublereal critTemperature() const;
/// critical pressure
virtual doublereal critPressure() const;
/// critical pressure
virtual doublereal critPressure() const;
/// critical density
virtual doublereal critDensity() const;
/// critical density
virtual doublereal critDensity() const;
/// saturation temperature
//virtual doublereal satTemperature(doublereal p) const;
/// saturation temperature
//virtual doublereal satTemperature(doublereal p) const;
/// saturation pressure
virtual doublereal satPressure(doublereal t);
/// saturation pressure
virtual doublereal satPressure(doublereal t);
virtual void setDensity(double dens);
double density() const;
virtual void setTemperature(double temp);
double temperature() const;
virtual void setState_TP(double temp, double pres);
virtual void setDensity(double dens);
double density() const;
virtual void setTemperature(double temp);
double temperature() const;
virtual void setState_TP(double temp, double pres);
doublereal molecularWeight() const;
void setMolecularWeight(double mw);
doublereal molecularWeight() const;
void setMolecularWeight(double mw);
virtual void constructPDSS(ThermoPhase *tp, int spindex);
virtual void constructPDSSFile(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id);
virtual void constructPDSSXML(ThermoPhase *tp, int spindex,
XML_Node& phaseNode, std::string id);
virtual void initThermoXML(XML_Node& eosdata, std::string id);
virtual void initThermo();
virtual void setParametersFromXML(const XML_Node& eosdata);
WaterPropsIAPWS *getWater() const {
return m_sub;
}
protected:
virtual void constructPDSS(ThermoPhase *tp, int spindex);
virtual void constructPDSSFile(ThermoPhase *tp, int spindex,
std::string inputFile, std::string id);
virtual void constructPDSSXML(ThermoPhase *tp, int spindex,
XML_Node& phaseNode, std::string id);
virtual void initThermoXML(XML_Node& eosdata, std::string id);
virtual void initThermo();
virtual void setParametersFromXML(const XML_Node& eosdata);
WaterPropsIAPWS *getWater() const {
return m_sub;
}
protected:
private:
mutable WaterPropsIAPWS *m_sub;
private:
mutable WaterPropsIAPWS *m_sub;
/**
* state of the system (temperature and density);
*/
doublereal m_temp;
doublereal m_dens;
/**
* state of the system (temperature and density);
*/
doublereal m_temp;
doublereal m_dens;
/*
* state of the fluid
* 0 gas
* 1 liquid
* 2 supercrit
*/
int m_iState;
/*
* state of the fluid
* 0 gas
* 1 liquid
* 2 supercrit
*/
int m_iState;
/**
* Offset constants used to obtain consistency with the NIST database.
* This is added to all internal energy and enthalpy results.
* units = J kmol-1.
*/
double EW_Offset;
/**
* Offset constants used to obtain consistency with the NIST database.
* This is added to all internal energy and enthalpy results.
* units = J kmol-1.
*/
double EW_Offset;
/*
* Offset constant used to obtain consistency with NIST convention.
* This is added to all internal entropy results.
* units = J kmol-1 K-1.
*/
double SW_Offset;
/*
* Offset constant used to obtain consistency with NIST convention.
* This is added to all internal entropy results.
* units = J kmol-1 K-1.
*/
double SW_Offset;
bool m_verbose;
bool m_verbose;
/**
* Since this phase represents a liquid phase, it's an error to
* return a gas-phase answer. However, if the below is true, then
* a gas-phase answer is allowed. This is used to check the thermodynamic
* consistency with ideal-gas thermo functions for example.
*/
bool m_allowGasPhase;
};
/**
* Since this phase represents a liquid phase, it's an error to
* return a gas-phase answer. However, if the below is true, then
* a gas-phase answer is allowed. This is used to check the thermodynamic
* consistency with ideal-gas thermo functions for example.
*/
bool m_allowGasPhase;
};
}

View file

@ -499,8 +499,21 @@ namespace Cantera {
}
setDensity(dd);
}
// Return the volumetric thermal expansion coefficient. Units: 1/K.
/*
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
doublereal WaterSSTP::thermalExpansionCoeff() const {
doublereal pres = pressure();
double T = temperature();
doublereal val = m_sub->coeffThermExp(T, pres);
return val;
}
// critical temperature
doublereal WaterSSTP::critTemperature() const { return m_sub->Tcrit(); }

View file

@ -186,6 +186,16 @@ namespace Cantera {
virtual doublereal pressure() const;
virtual void setPressure(doublereal p);
//! Return the volumetric thermal expansion coefficient. Units: 1/K.
/*!
* The thermal expansion coefficient is defined as
* \f[
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
* \f]
*/
virtual doublereal thermalExpansionCoeff() const;
/**
* @}
* @name Potential Energy