Merged changes from the trunk.
This commit is contained in:
Harry Moffat
parent
1687ecaf1b
commit
6b1a1ac27f
10 changed files with 1450 additions and 317 deletions
|
|
@ -252,7 +252,7 @@ namespace Cantera {
|
|||
*
|
||||
* where we can use the concept of microscopic reversibility to
|
||||
* write the reverse rate constant in terms of the
|
||||
* forward reate constant and the concentration equilibrium
|
||||
* forward rate constant and the concentration equilibrium
|
||||
* constant, \f$ K_c \f$.
|
||||
*
|
||||
* \f[
|
||||
|
|
|
|||
|
|
@ -66,7 +66,7 @@ ELECTRO_OBJ = MolalityVPSSTP.o VPStandardStateTP.o \
|
|||
WaterPropsIAPWSphi.o WaterPropsIAPWS.o WaterProps.o \
|
||||
PDSS.o PDSS_Water.o PDSS_HKFT.o \
|
||||
HMWSoln.o HMWSoln_input.o DebyeHuckel.o \
|
||||
WaterSSTP.o \
|
||||
WaterSSTP.o MetalSHEelectrons.o \
|
||||
VPSSMgr_Water_ConstVol.o VPSSMgr_Water_HKFT.o
|
||||
|
||||
ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \
|
||||
|
|
@ -74,7 +74,7 @@ ELECTRO_H = MolalityVPSSTP.h VPStandardStateTP.h \
|
|||
WaterPropsIAPWSphi.h WaterPropsIAPWS.h WaterProps.h \
|
||||
PDSS.h PDSS_Water.h PDSS_HKFT.h \
|
||||
HMWSoln.h electrolytes.h \
|
||||
DebyeHuckel.h WaterSSTP.h VPSSMgr_Water_HKFT.h \
|
||||
DebyeHuckel.h WaterSSTP.h MetalSHEelectrons.h VPSSMgr_Water_HKFT.h \
|
||||
VPSSMgr_Water_ConstVol.h
|
||||
endif
|
||||
ifeq ($(do_issp),1)
|
||||
|
|
|
|||
546
Cantera/src/thermo/MetalSHEelectrons.cpp
Normal file
546
Cantera/src/thermo/MetalSHEelectrons.cpp
Normal file
|
|
@ -0,0 +1,546 @@
|
|||
/**
|
||||
* @file MetalSHEelectrons.cpp
|
||||
* Definition file for the %MetalSHEElectrons class, which represents the
|
||||
* electrons in a metal that are consistent with the
|
||||
* SHE electrode (see \ref thermoprops and
|
||||
* class \link Cantera::MetalSHEelectrons MetalSHEelectrons\endlink)
|
||||
*/
|
||||
|
||||
/*
|
||||
* Copywrite (2005) Sandia Corporation. Under the terms of
|
||||
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
|
||||
* U.S. Government retains certain rights in this software.
|
||||
*
|
||||
*/
|
||||
|
||||
/*
|
||||
* $Id: MetalSHEelectrons.cpp 279 2009-12-05 19:08:43Z hkmoffa $
|
||||
*/
|
||||
|
||||
#include "ct_defs.h"
|
||||
|
||||
#include "MetalSHEelectrons.h"
|
||||
#include "SingleSpeciesTP.h"
|
||||
#include "ThermoFactory.h"
|
||||
|
||||
#include <string>
|
||||
using namespace std;
|
||||
using namespace Cantera;
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/*
|
||||
* ---- Constructors -------
|
||||
*/
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Default Constructor for the MetalSHEelectrons class
|
||||
*/
|
||||
MetalSHEelectrons::MetalSHEelectrons():
|
||||
SingleSpeciesTP(),
|
||||
xdef_(0)
|
||||
{
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Create and initialize a MetalSHEelectrons ThermoPhase object
|
||||
// from an asci input file
|
||||
/*
|
||||
* @param infile name of the input file
|
||||
* @param id name of the phase id in the file.
|
||||
* If this is blank, the first phase in the file is used.
|
||||
*/
|
||||
MetalSHEelectrons::MetalSHEelectrons(std::string infile, std::string id) :
|
||||
SingleSpeciesTP(),
|
||||
xdef_(0)
|
||||
{
|
||||
XML_Node* root;
|
||||
if (infile == "MetalSHEelectrons_default.xml") {
|
||||
xdef_ = MetalSHEelectrons::makeDefaultXMLTree();
|
||||
root = xdef_;
|
||||
} else {
|
||||
root = get_XML_File(infile);
|
||||
}
|
||||
if (id == "-") id = "";
|
||||
XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root);
|
||||
if (!xphase) {
|
||||
throw CanteraError("MetalSHEelectrons::MetalSHEelectrons",
|
||||
"Couldn't find phase name in file:" + id);
|
||||
}
|
||||
// Check the model name to ensure we have compatibility
|
||||
const XML_Node& th = xphase->child("thermo");
|
||||
std::string model = th["model"];
|
||||
if (model != "MetalSHEelectrons") {
|
||||
throw CanteraError("MetalSHEelectrons::MetalSHEelectrons",
|
||||
"thermo model attribute must be MetalSHEelectrons");
|
||||
}
|
||||
importPhase(*xphase, this);
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Full Constructor.
|
||||
/*
|
||||
* @param phaseRef XML node pointing to a MetalSHEelectrons description
|
||||
* @param id Id of the phase.
|
||||
*/
|
||||
MetalSHEelectrons::MetalSHEelectrons(XML_Node& xmlphase, std::string id) :
|
||||
SingleSpeciesTP(),
|
||||
xdef_(0)
|
||||
{
|
||||
if (id != "") {
|
||||
std::string idxml = xmlphase["id"];
|
||||
if (id != idxml) {
|
||||
throw CanteraError("MetalSHEelectrons::MetalSHEelectrons",
|
||||
"id's don't match");
|
||||
}
|
||||
}
|
||||
const XML_Node& th = xmlphase.child("thermo");
|
||||
std::string model = th["model"];
|
||||
if (model != "MetalSHEelectrons") {
|
||||
throw CanteraError("MetalSHEelectrons::MetalSHEelectrons",
|
||||
"thermo model attribute must be MetalSHEelectrons");
|
||||
}
|
||||
importPhase(xmlphase, this);
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Copy constructor
|
||||
/*
|
||||
* @param right Object to be copied
|
||||
*/
|
||||
MetalSHEelectrons::MetalSHEelectrons(const MetalSHEelectrons &right) :
|
||||
SingleSpeciesTP()
|
||||
{
|
||||
operator=(right);
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Destructor for the routine (virtual)
|
||||
*
|
||||
*/
|
||||
MetalSHEelectrons::~MetalSHEelectrons()
|
||||
{
|
||||
if (xdef_) {
|
||||
delete xdef_;
|
||||
}
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Assignment operator
|
||||
/*
|
||||
* @param right Object to be copied
|
||||
*/
|
||||
MetalSHEelectrons &
|
||||
MetalSHEelectrons::operator=(const MetalSHEelectrons & right)
|
||||
{
|
||||
if (&right != this) {
|
||||
SingleSpeciesTP::operator=(right);
|
||||
}
|
||||
|
||||
if (xdef_) {
|
||||
delete xdef_;
|
||||
}
|
||||
xdef_ = new XML_Node(*right.xdef_);
|
||||
|
||||
return *this;
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Duplication function
|
||||
/*
|
||||
* This virtual function is used to create a duplicate of the
|
||||
* current phase. It's used to duplicate the phase when given
|
||||
* a ThermoPhase pointer to the phase.
|
||||
*
|
||||
* @return It returns a ThermoPhase pointer.
|
||||
*/
|
||||
ThermoPhase *MetalSHEelectrons::duplMyselfAsThermoPhase() const {
|
||||
MetalSHEelectrons *stp = new MetalSHEelectrons(*this);
|
||||
return (ThermoPhase *) stp;
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
/*
|
||||
* ---- Utilities -----
|
||||
*/
|
||||
|
||||
/*
|
||||
* Equation of state flag. Returns the value cStoichSubstance,
|
||||
* defined in mix_defs.h.
|
||||
*/
|
||||
int MetalSHEelectrons::eosType() const {
|
||||
return cMetalSHEelectrons;
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
/*
|
||||
* ---- Molar Thermodynamic properties of the solution ----
|
||||
*/
|
||||
|
||||
/**
|
||||
* ----- Mechanical Equation of State ------
|
||||
*/
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Pressure. Units: Pa.
|
||||
* For an incompressible substance, the density is independent
|
||||
* of pressure. This method simply returns the stored
|
||||
* pressure value.
|
||||
*/
|
||||
doublereal MetalSHEelectrons::pressure() const {
|
||||
return m_press;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Set the pressure at constant temperature. Units: Pa.
|
||||
* For an incompressible substance, the density is
|
||||
* independent of pressure. Therefore, this method only
|
||||
* stores the specified pressure value. It does not
|
||||
* modify the density.
|
||||
*/
|
||||
void MetalSHEelectrons::setPressure(doublereal p) {
|
||||
m_press = p;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* The isothermal compressibility. Units: 1/Pa.
|
||||
* The isothermal compressibility is defined as
|
||||
* \f[
|
||||
* \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
|
||||
* \f]
|
||||
*
|
||||
* It's equal to zero for this model, since the molar volume
|
||||
* doesn't change with pressure or temperature.
|
||||
*/
|
||||
doublereal MetalSHEelectrons::isothermalCompressibility() const {
|
||||
return -1.0/pressure();
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* The 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]
|
||||
*
|
||||
* It's equal to zero for this model, since the molar volume
|
||||
* doesn't change with pressure or temperature.
|
||||
*/
|
||||
doublereal MetalSHEelectrons::thermalExpansionCoeff() const {
|
||||
return 1.0/temperature();
|
||||
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* ---- Chemical Potentials and Activities ----
|
||||
*/
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* This method returns the array of generalized
|
||||
* concentrations. For a stoichiomeetric substance, there is
|
||||
* only one species, and the generalized concentration is 1.0.
|
||||
*/
|
||||
void MetalSHEelectrons::
|
||||
getActivityConcentrations(doublereal* c) const {
|
||||
c[0] = 1.0;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* The standard concentration. This is defined as the concentration
|
||||
* by which the generalized concentration is normalized to produce
|
||||
* the activity.
|
||||
*/
|
||||
doublereal MetalSHEelectrons::standardConcentration(int k) const {
|
||||
return 1.0;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Returns the natural logarithm of the standard
|
||||
* concentration of the kth species
|
||||
*/
|
||||
doublereal MetalSHEelectrons::logStandardConc(int k) const {
|
||||
return 0.0;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Returns the units of the standard and generalized
|
||||
* concentrations Note they have the same units, as their
|
||||
* ratio is defined to be equal to the activity of the kth
|
||||
* species in the solution, which is unitless.
|
||||
*
|
||||
* This routine is used in print out applications where the
|
||||
* units are needed. Usually, MKS units are assumed throughout
|
||||
* the program and in the XML input files.
|
||||
*
|
||||
* uA[0] = kmol units - default = 1
|
||||
* uA[1] = m units - default = -nDim(), the number of spatial
|
||||
* dimensions in the Phase class.
|
||||
* uA[2] = kg units - default = 0;
|
||||
* uA[3] = Pa(pressure) units - default = 0;
|
||||
* uA[4] = Temperature units - default = 0;
|
||||
* uA[5] = time units - default = 0
|
||||
*/
|
||||
void MetalSHEelectrons::
|
||||
getUnitsStandardConc(doublereal *uA, int k, int sizeUA) const {
|
||||
for (int i = 0; i < 6; i++) {
|
||||
uA[i] = 0;
|
||||
}
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* ---- Partial Molar Properties of the Solution ----
|
||||
*/
|
||||
|
||||
//====================================================================================================================
|
||||
|
||||
/*
|
||||
* ---- Properties of the Standard State of the Species in the Solution
|
||||
* ----
|
||||
*/
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Get the array of chemical potentials at unit activity
|
||||
* \f$ \mu^0_k \f$.
|
||||
*
|
||||
* For a stoichiometric substance, there is no activity term in
|
||||
* the chemical potential expression, and therefore the
|
||||
* standard chemical potential and the chemical potential
|
||||
* are both equal to the molar Gibbs function.
|
||||
*/
|
||||
void MetalSHEelectrons::
|
||||
getStandardChemPotentials(doublereal* mu0) const {
|
||||
getGibbs_RT(mu0);
|
||||
mu0[0] *= GasConstant * temperature();
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Get the nondimensional Enthalpy functions for the species
|
||||
* at their standard states at the current
|
||||
* <I>T</I> and <I>P</I> of the solution.
|
||||
* Molar enthalpy. Units: J/kmol. For an incompressible,
|
||||
* stoichiometric substance, the internal energy is
|
||||
* independent of pressure, and therefore the molar enthalpy
|
||||
* is \f[ \hat h(T, P) = \hat u(T) + P \hat v \f], where the
|
||||
* molar specific volume is constant.
|
||||
*/
|
||||
void MetalSHEelectrons::getEnthalpy_RT(doublereal* hrt) const {
|
||||
getEnthalpy_RT_ref(hrt);
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Get the array of nondimensional Entropy functions for the
|
||||
* standard state species
|
||||
* at the current <I>T</I> and <I>P</I> of the solution.
|
||||
*/
|
||||
void MetalSHEelectrons::getEntropy_R(doublereal* sr) const {
|
||||
getEntropy_R_ref(sr);
|
||||
doublereal tmp = log (pressure() / m_p0);
|
||||
sr[0] -= tmp;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Get the nondimensional Gibbs functions for the species
|
||||
* at their standard states of solution at the current T and P
|
||||
* of the solution
|
||||
*/
|
||||
void MetalSHEelectrons::getGibbs_RT(doublereal* grt) const {
|
||||
getGibbs_RT_ref(grt);
|
||||
doublereal tmp = log (pressure() / m_p0);
|
||||
grt[0] += tmp;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Get the nondimensional Gibbs functions for the standard
|
||||
* state of the species at the current T and P.
|
||||
*/
|
||||
void MetalSHEelectrons::getCp_R(doublereal* cpr) const {
|
||||
_updateThermo();
|
||||
cpr[0] = m_cp0_R[0];
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Molar internal energy (J/kmol).
|
||||
* For an incompressible,
|
||||
* stoichiometric substance, the molar internal energy is
|
||||
* independent of pressure. Since the thermodynamic properties
|
||||
* are specified by giving the standard-state enthalpy, the
|
||||
* term \f$ P_0 \hat v\f$ is subtracted from the specified molar
|
||||
* enthalpy to compute the molar internal energy.
|
||||
*/
|
||||
void MetalSHEelectrons::getIntEnergy_RT(doublereal* urt) const {
|
||||
getEnthalpy_RT(urt);
|
||||
urt[0] -= 1.0;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* ---- Thermodynamic Values for the Species Reference States ----
|
||||
*/
|
||||
/*
|
||||
* Molar internal energy or the reference state at the current
|
||||
* temperature, T (J/kmol).
|
||||
* For an incompressible,
|
||||
* stoichiometric substance, the molar internal energy is
|
||||
* independent of pressure. Since the thermodynamic properties
|
||||
* are specified by giving the standard-state enthalpy, the
|
||||
* term \f$ P_0 \hat v\f$ is subtracted from the specified molar
|
||||
* enthalpy to compute the molar internal energy.
|
||||
*
|
||||
* Note, this is equal to the standard state internal energy
|
||||
* evaluated at the reference pressure.
|
||||
*/
|
||||
void MetalSHEelectrons::getIntEnergy_RT_ref(doublereal* urt) const {
|
||||
_updateThermo();
|
||||
doublereal RT = GasConstant * temperature();
|
||||
doublereal PV = m_p0 / molarDensity();
|
||||
urt[0] = m_h0_RT[0] - PV / RT;
|
||||
}
|
||||
|
||||
/*
|
||||
* ---- Saturation Properties
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* ---- Initialization and Internal functions
|
||||
*/
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* @internal Initialize. This method is provided to allow
|
||||
* subclasses to perform any initialization required after all
|
||||
* species have been added. For example, it might be used to
|
||||
* resize internal work arrays that must have an entry for
|
||||
* each species. The base class implementation does nothing,
|
||||
* and subclasses that do not require initialization do not
|
||||
* need to overload this method. When importing a CTML phase
|
||||
* description, this method is called just prior to returning
|
||||
* from function importPhase.
|
||||
*
|
||||
* @see importCTML.cpp
|
||||
*/
|
||||
void MetalSHEelectrons::initThermo() {
|
||||
/*
|
||||
* Call the base class thermo initializer
|
||||
*/
|
||||
SingleSpeciesTP::initThermo();
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
void MetalSHEelectrons::initThermoXML(XML_Node& phaseNode, std::string id) {
|
||||
/*
|
||||
* Find the Thermo XML node
|
||||
*/
|
||||
if (!phaseNode.hasChild("thermo")) {
|
||||
throw CanteraError("MetalSHEelectrons::initThermoXML",
|
||||
"no thermo XML node");
|
||||
}
|
||||
XML_Node &tnode = phaseNode.child("thermo");
|
||||
doublereal dens = 2.65E3;
|
||||
if (tnode.hasChild("density")) {
|
||||
dens = getFloatDefaultUnits(tnode, "density", "kg/m3");
|
||||
}
|
||||
setDensity(dens);
|
||||
SingleSpeciesTP::initThermoXML(phaseNode, id);
|
||||
}
|
||||
//====================================================================================================================
|
||||
XML_Node *MetalSHEelectrons::makeDefaultXMLTree()
|
||||
{
|
||||
XML_Node *xtop = new XML_Node("ctml", 0);
|
||||
XML_Node &xv = xtop->addChild("validate");
|
||||
xv.addAttribute("reactions", "yes");
|
||||
xv.addAttribute("species", "yes");
|
||||
|
||||
XML_Node &xp = xtop->addChild("phase");
|
||||
xp.addAttribute("dim", "3");
|
||||
xp.addAttribute("id", "MetalSHEelectrons");
|
||||
XML_Node &xe = xp.addChild("elementArray", "E");
|
||||
xe.addAttribute("datasrc", "elements.xml");
|
||||
XML_Node &xs = xp.addChild("speciesArray", "she_electron");
|
||||
xs.addAttribute("datasrc", "#species_Metal_SHEelectrons");
|
||||
XML_Node &xt = xp.addChild("thermo");
|
||||
xt.addAttribute("model", "metalSHEelectrons");
|
||||
XML_Node &xtr = xp.addChild("transport");
|
||||
xtr.addAttribute("model", "none");
|
||||
XML_Node &xk = xp.addChild("kinetics");
|
||||
xk.addAttribute("model", "none");
|
||||
|
||||
XML_Node &xsd = xtop->addChild("speciesData");
|
||||
xsd.addAttribute("id", "species_Metal_SHEelectrons");
|
||||
|
||||
XML_Node &xsp = xsd.addChild("species");
|
||||
xsp.addAttribute("name", "she_electron");
|
||||
xsp.addChild("atomArray", "E:1");
|
||||
xsp.addChild("charge", "-1");
|
||||
XML_Node &xspt = xsp.addChild("thermo");
|
||||
|
||||
XML_Node &xN1 = xspt.addChild("NASA");
|
||||
xN1.addAttribute("Tmax", "1000.");
|
||||
xN1.addAttribute("Tmin", "200.");
|
||||
xN1.addAttribute("P0", "100000.0");
|
||||
XML_Node &xF1 = xsd.addChild("floatArray",
|
||||
"1.172165560E+00, 3.990260375E-03, -9.739075500E-06, "
|
||||
"1.007860470E-08, -3.688058805E-12, -4.589675865E+02, 3.415051190E-01" );
|
||||
xF1.addAttribute("name", "coeffs");
|
||||
xF1.addAttribute("size", "7");
|
||||
|
||||
XML_Node &xN2 = xspt.addChild("NASA");
|
||||
xN2.addAttribute("Tmax", "6000.");
|
||||
xN2.addAttribute("Tmin", "1000.");
|
||||
xN2.addAttribute("P0", "100000.0");
|
||||
XML_Node &xF2 = xsd.addChild("floatArray",
|
||||
"1.466432895E+00, 4.133039835E-04, -7.320116750E-08, 7.705017950E-12,"
|
||||
"-3.444022160E-16, -4.065327985E+02, -5.121644350E-01");
|
||||
xF2.addAttribute("name", "coeffs");
|
||||
xF2.addAttribute("size", "7");
|
||||
|
||||
return xtop;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* setParameters:
|
||||
*
|
||||
* Generic routine that is used to set the parameters used
|
||||
* by this model.
|
||||
* C[0] = density of phase [ kg/m3 ]
|
||||
*/
|
||||
void MetalSHEelectrons::setParameters(int n, doublereal * const c) {
|
||||
doublereal rho = c[0];
|
||||
setDensity(rho);
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* getParameters:
|
||||
*
|
||||
* Generic routine that is used to get the parameters used
|
||||
* by this model.
|
||||
* n = 1
|
||||
* C[0] = density of phase [ kg/m3 ]
|
||||
*/
|
||||
void MetalSHEelectrons::getParameters(int &n, doublereal * const c) const {
|
||||
doublereal rho = density();
|
||||
n = 1;
|
||||
c[0] = rho;
|
||||
}
|
||||
//====================================================================================================================
|
||||
/*
|
||||
* Reads an xml data block for the parameters needed by this
|
||||
* routine. eosdata is a reference to the xml thermo block, and looks
|
||||
* like this:
|
||||
*
|
||||
* <phase id="stoichsolid" >
|
||||
* <thermo model="StoichSubstance">
|
||||
* <density units="g/cm3">3.52</density>
|
||||
* </thermo>
|
||||
* </phase>
|
||||
*/
|
||||
void MetalSHEelectrons::setParametersFromXML(const XML_Node& eosdata) {
|
||||
std::string model = eosdata["model"];
|
||||
if (model != "MetalSHEelectrons") {
|
||||
throw CanteraError("MetalSHEelectrons::setParametersFromXML",
|
||||
"thermo model attribute must be MetalSHEelectrons");
|
||||
}
|
||||
doublereal rho = 2.65E3;
|
||||
if (eosdata.hasChild("density")) {
|
||||
rho = getFloat(eosdata, "density", "toSI");
|
||||
}
|
||||
setDensity(rho);
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
}
|
||||
565
Cantera/src/thermo/MetalSHEelectrons.h
Normal file
565
Cantera/src/thermo/MetalSHEelectrons.h
Normal file
|
|
@ -0,0 +1,565 @@
|
|||
/**
|
||||
* @file MetalSHEelectrons.h
|
||||
* Header file for the %MetalSHEElectrons class, which represents the
|
||||
* electrons in a metal that are consistent with the
|
||||
* SHE electrode (see \ref thermoprops and
|
||||
* class \link Cantera::MetalSHEelectrons MetalSHEelectrons\endlink)
|
||||
*/
|
||||
|
||||
/*
|
||||
* Copywrite (2005) Sandia Corporation. Under the terms of
|
||||
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
|
||||
* U.S. Government retains certain rights in this software.
|
||||
*/
|
||||
|
||||
/*
|
||||
* $Date: 2009-12-05 12:08:43 -0700 (Sat, 05 Dec 2009) $
|
||||
* $Revision: 279 $
|
||||
*/
|
||||
|
||||
#ifndef CT_METALSHEELECTRONS_H
|
||||
#define CT_METALSHEELECTRONS_H
|
||||
|
||||
#include "mix_defs.h"
|
||||
#include "SingleSpeciesTP.h"
|
||||
//#include "SpeciesThermo.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
//! Class %MetalSHEelectrons represents electrons within
|
||||
//! a metal, adjacent to an aqueous electrolyte, that are consistent with the SHE reference electrode.
|
||||
/*!
|
||||
* The class is based on the electron having a chemical potential
|
||||
* equal to one-half of the entropy of the H<SUP>2</SUP> gas at the system pressure
|
||||
*
|
||||
*
|
||||
* <b> Specification of Species Standard %State Properties </b>
|
||||
*
|
||||
* This class inherits from SingleSpeciesTP.
|
||||
* It is assumed that the reference state thermodynamics may be
|
||||
* obtained by a pointer to a populated species thermodynamic property
|
||||
* manager class (see ThermoPhase::m_spthermo). How to relate pressure
|
||||
* changes to the reference state thermodynamics is resolved at this level.
|
||||
*
|
||||
* The enthalpy function is given by the following relation.
|
||||
*
|
||||
* \f[
|
||||
* h^o_k(T,P) = h^{ref}_k(T)
|
||||
* \f]
|
||||
*
|
||||
* The standard state constant-pressure heat capacity is independent of pressure:
|
||||
*
|
||||
* \f[
|
||||
* Cp^o_k(T,P) = Cp^{ref}_k(T)
|
||||
* \f]
|
||||
*
|
||||
* The standard state entropy depends in the following fashion on pressure:
|
||||
*
|
||||
* \f[
|
||||
* S^o_k(T,P) = S^{ref}_k(T) - R \ln(\frac{P}{P_{ref}})
|
||||
* \f]
|
||||
*
|
||||
* The standard state gibbs free energy is obtained from the enthalpy and entropy
|
||||
* functions:
|
||||
*
|
||||
* \f[
|
||||
* \mu^o_k(T,P) = h^o_k(T,P) - S^o_k(T,P) T
|
||||
* \f]
|
||||
*
|
||||
* \f[
|
||||
* \mu^o_k(T,P) = \mu^{ref}_k(T) + R T \ln( \frac{P}{P_{ref}})
|
||||
* \f]
|
||||
*
|
||||
* where
|
||||
* \f[
|
||||
* \mu^{ref}_k(T) = h^{ref}_k(T) - T S^{ref}_k(T)
|
||||
* \f]
|
||||
*
|
||||
* The standard state internal energy is obtained from the enthalpy function also
|
||||
*
|
||||
* \f[
|
||||
* u^o_k(T,P) = h^o_k(T) - R T
|
||||
* \f]
|
||||
*
|
||||
*
|
||||
* <b> Specification of Solution Thermodynamic Properties </b>
|
||||
*
|
||||
* All solution properties are obtained from the standard state
|
||||
* species functions, since there is only one species in the phase.
|
||||
*
|
||||
* <b> %Application within %Kinetics Managers </b>
|
||||
*
|
||||
* The standard concentration is equal to 1.0. This means that the
|
||||
* kinetics operator works on an activities basis. Since this
|
||||
* is a stoichiometric substance, this means that the concentration
|
||||
* of this phase drops out of kinetics expressions since the activity is
|
||||
* always equal to one.
|
||||
*
|
||||
* This is what is expected of electrons. The only effect that this class will
|
||||
* have on reactions is in terms of the standard state chemical potential, which
|
||||
* is equal to 1/2 of the H2 gas chemical potential, and the voltage assigned
|
||||
* to the electron, which is the voltage of the metal.
|
||||
*
|
||||
*
|
||||
* <b> Instanteation of the Class </b>
|
||||
*
|
||||
* The constructor for this phase is located in the default ThermoFactory
|
||||
* for %Cantera. A new %MetalSHEelectrons object may be created by
|
||||
* the following code snippets, where the file metalSHEelectrons.xml exists
|
||||
* in a local directory:
|
||||
*
|
||||
* @code
|
||||
* MetalSHEelectrons *eMetal = new MetalSHEelectrons("metalSHEelectrons.xml", "");
|
||||
* @endcode
|
||||
*
|
||||
* or by the following call to importPhase():
|
||||
*
|
||||
* @code
|
||||
* sprintf(file_ID,"%s#MetalSHEelectrons", iFile);
|
||||
* XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
|
||||
* MetalSHEelectrons eMetal;
|
||||
* importPhase(*xm, &eMetal);
|
||||
* @endcode
|
||||
*
|
||||
* @code
|
||||
* ThermoPhase *eMetal = newPhase(" MetalSHEelectrons.xml", "MetalSHEelectrons");
|
||||
* @endcode
|
||||
*
|
||||
* Additionally, this phase may be created without including an xml file with
|
||||
* the special command, where the default file is embedded into this object.
|
||||
*
|
||||
* @code
|
||||
* MetalSHEelectrons *eMetal = new MetalSHEelectrons("MetalSHEelectrons_default.xml", "");
|
||||
* @endcode
|
||||
*
|
||||
*
|
||||
*
|
||||
* <b> XML Example </b>
|
||||
*
|
||||
* The phase model name for this is called %MetalSHEelectrons. It must be supplied
|
||||
* as the model attribute of the thermo XML element entry.
|
||||
* Within the phase XML block,
|
||||
* the density of the phase must be specified though it's not used. An example of an XML file
|
||||
* this phase is given below.
|
||||
*
|
||||
* @verbatim
|
||||
<?xml version="1.0"?>
|
||||
<ctml>
|
||||
<validate reactions="yes" species="yes"/>
|
||||
|
||||
<phase dim="3" id="MetalSHEelectrons">
|
||||
<elementArray datasrc="elements.xml">
|
||||
E
|
||||
</elementArray>
|
||||
<speciesArray datasrc="#species_Metal_SHEelectrons"> she_electron </speciesArray>
|
||||
<thermo model="metalSHEelectrons">
|
||||
<density units="g/cm3">2.165</density>
|
||||
</thermo>
|
||||
<transport model="None"/>
|
||||
<kinetics model="none"/>
|
||||
</phase>
|
||||
|
||||
<!-- species definitions -->
|
||||
<speciesData id="species_Metal_SHEelectrons">
|
||||
<species name="she_electron">
|
||||
<atomArray> E:1 </atomArray>
|
||||
<charge> -1 </charge>
|
||||
<thermo>
|
||||
<NASA Tmax="1000.0" Tmin="200.0" P0="100000.0">
|
||||
<floatArray name="coeffs" size="7">
|
||||
1.172165560E+00, 3.990260375E-03, -9.739075500E-06, 1.007860470E-08,
|
||||
-3.688058805E-12, -4.589675865E+02, 3.415051190E-01
|
||||
</floatArray>
|
||||
</NASA>
|
||||
<NASA Tmax="6000.0" Tmin="1000.0" P0="100000.0">
|
||||
<floatArray name="coeffs" size="7">
|
||||
1.466432895E+00, 4.133039835E-04, -7.320116750E-08, 7.705017950E-12,
|
||||
-3.444022160E-16, -4.065327985E+02, -5.121644350E-01
|
||||
</floatArray>
|
||||
</NASA>
|
||||
</thermo>
|
||||
<density units="g/cm3">2.165</density>
|
||||
</species>
|
||||
</speciesData>
|
||||
</ctml>
|
||||
@endverbatim
|
||||
*
|
||||
* The model attribute, "MetalSHEelectrons", on the thermo element
|
||||
* identifies the phase as being a %MetalSHEelectrons object.
|
||||
*
|
||||
* @ingroup thermoprops
|
||||
*/
|
||||
class MetalSHEelectrons : public SingleSpeciesTP {
|
||||
|
||||
public:
|
||||
|
||||
//! Default constructor for the MetalSHEelectrons class
|
||||
MetalSHEelectrons();
|
||||
|
||||
//! Construct and initialize a %MetalSHEelectrons %ThermoPhase object
|
||||
//! directly from an asci input file
|
||||
/*!
|
||||
* @param infile name of the input file
|
||||
* @param id name of the phase id in the file.
|
||||
* If this is blank, the first phase in the file is used.
|
||||
*/
|
||||
MetalSHEelectrons(std::string infile, std::string id = "");
|
||||
|
||||
//! Construct and initialize a MetalSHEelectrons ThermoPhase object
|
||||
//! directly from an XML database
|
||||
/*!
|
||||
* @param phaseRef XML node pointing to a MetalSHEelectrons description
|
||||
* @param id Id of the phase.
|
||||
*/
|
||||
MetalSHEelectrons(XML_Node& phaseRef, std::string id = "");
|
||||
|
||||
//! Copy constructor
|
||||
/*!
|
||||
* @param right Object to be copied
|
||||
*/
|
||||
MetalSHEelectrons(const MetalSHEelectrons &right);
|
||||
|
||||
//! Assignment operator
|
||||
/*!
|
||||
* @param right Object to be copied
|
||||
*/
|
||||
MetalSHEelectrons & operator=(const MetalSHEelectrons & right);
|
||||
|
||||
//! Destructor for the routine (virtual)
|
||||
virtual ~MetalSHEelectrons();
|
||||
|
||||
//! Duplication function
|
||||
/*!
|
||||
* This virtual function is used to create a duplicate of the
|
||||
* current phase. It's used to duplicate the phase when given
|
||||
* a ThermoPhase pointer to the phase.
|
||||
*
|
||||
* @return It returns a ThermoPhase pointer.
|
||||
*/
|
||||
ThermoPhase *duplMyselfAsThermoPhase() const;
|
||||
|
||||
/**
|
||||
*
|
||||
* @name Utilities
|
||||
* @{
|
||||
*/
|
||||
|
||||
/**
|
||||
* Equation of state flag.
|
||||
*
|
||||
* Returns the value cStoichSubstance, defined in mix_defs.h.
|
||||
*/
|
||||
virtual int eosType() const;
|
||||
|
||||
/**
|
||||
* @}
|
||||
* @name Molar Thermodynamic Properties of the Solution
|
||||
* @{
|
||||
*/
|
||||
|
||||
/**
|
||||
* @}
|
||||
* @name Mechanical Equation of State
|
||||
* @{
|
||||
*/
|
||||
|
||||
|
||||
//! Report the Pressure. Units: Pa.
|
||||
/*!
|
||||
* For an incompressible substance, the density is independent
|
||||
* of pressure. This method simply returns the storred
|
||||
* pressure value.
|
||||
*/
|
||||
virtual doublereal pressure() const;
|
||||
|
||||
//! Set the pressure at constant temperature. Units: Pa.
|
||||
/*!
|
||||
* For an incompressible substance, the density is
|
||||
* independent of pressure. Therefore, this method only
|
||||
* stores the specified pressure value. It does not
|
||||
* modify the density.
|
||||
*
|
||||
* @param p Pressure (units - Pa)
|
||||
*/
|
||||
virtual void setPressure(doublereal p);
|
||||
|
||||
//! Returns the isothermal compressibility. Units: 1/Pa.
|
||||
/*!
|
||||
* The isothermal compressibility is defined as
|
||||
* \f[
|
||||
* \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
|
||||
* \f]
|
||||
*/
|
||||
virtual doublereal isothermalCompressibility() const;
|
||||
|
||||
//! 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 Activities, Standard States, and Activity Concentrations
|
||||
*
|
||||
* This section is largely handled by parent classes, since there
|
||||
* is only one species. Therefore, the activity is equal to one.
|
||||
* @{
|
||||
*/
|
||||
|
||||
//! This method returns an array of generalized concentrations
|
||||
/*!
|
||||
* \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k /
|
||||
* C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration
|
||||
* defined below and \f$ a_k \f$ are activities used in the
|
||||
* thermodynamic functions. These activity (or generalized)
|
||||
* concentrations are used
|
||||
* by kinetics manager classes to compute the forward and
|
||||
* reverse rates of elementary reactions.
|
||||
*
|
||||
* For a stoichiomeetric substance, there is
|
||||
* only one species, and the generalized concentration is 1.0.
|
||||
*
|
||||
* @param c Output array of generalized concentrations. The
|
||||
* units depend upon the implementation of the
|
||||
* reaction rate expressions within the phase.
|
||||
*/
|
||||
virtual void getActivityConcentrations(doublereal* c) const;
|
||||
|
||||
//! Return the standard concentration for the kth species
|
||||
/*!
|
||||
* The standard concentration \f$ C^0_k \f$ used to normalize
|
||||
* the activity (i.e., generalized) concentration.
|
||||
* This phase assumes that the kinetics operator works on an
|
||||
* dimensionless basis. Thus, the standard concentration is
|
||||
* equal to 1.0.
|
||||
*
|
||||
* @param k Optional parameter indicating the species. The default
|
||||
* is to assume this refers to species 0.
|
||||
* @return
|
||||
* Returns The standard Concentration as 1.0
|
||||
*/
|
||||
virtual doublereal standardConcentration(int k=0) const;
|
||||
|
||||
//! Natural logarithm of the standard concentration of the kth species.
|
||||
/*!
|
||||
* @param k index of the species (defaults to zero)
|
||||
*/
|
||||
virtual doublereal logStandardConc(int k=0) const;
|
||||
|
||||
//! Get the array of chemical potentials at unit activity for the species
|
||||
//! at their standard states at the current <I>T</I> and <I>P</I> of the solution.
|
||||
/*!
|
||||
* For a stoichiometric substance, there is no activity term in
|
||||
* the chemical potential expression, and therefore the
|
||||
* standard chemical potential and the chemical potential
|
||||
* are both equal to the molar Gibbs function.
|
||||
*
|
||||
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
|
||||
* \f$. The values are evaluated at the current
|
||||
* temperature and pressure of the solution
|
||||
*
|
||||
* @param mu0 Output vector of chemical potentials.
|
||||
* Length: m_kk.
|
||||
*/
|
||||
virtual void getStandardChemPotentials(doublereal* mu0) const;
|
||||
|
||||
//! Returns the units of the standard and generalized concentrations.
|
||||
/*!
|
||||
* Note they have the same units, as their
|
||||
* ratio is defined to be equal to the activity of the kth
|
||||
* species in the solution, which is unitless.
|
||||
*
|
||||
* This routine is used in print out applications where the
|
||||
* units are needed. Usually, MKS units are assumed throughout
|
||||
* the program and in the XML input files.
|
||||
*
|
||||
* The base %ThermoPhase class assigns thedefault quantities
|
||||
* of (kmol/m3) for all species.
|
||||
* Inherited classes are responsible for overriding the default
|
||||
* values if necessary.
|
||||
*
|
||||
* @param uA Output vector containing the units
|
||||
* uA[0] = kmol units - default = 1
|
||||
* uA[1] = m units - default = -nDim(), the number of spatial
|
||||
* dimensions in the Phase class.
|
||||
* uA[2] = kg units - default = 0;
|
||||
* uA[3] = Pa(pressure) units - default = 0;
|
||||
* uA[4] = Temperature units - default = 0;
|
||||
* uA[5] = time units - default = 0
|
||||
* @param k species index. Defaults to 0.
|
||||
* @param sizeUA output int containing the size of the vector.
|
||||
* Currently, this is equal to 6.
|
||||
*/
|
||||
virtual void getUnitsStandardConc(doublereal *uA, int k = 0,
|
||||
int sizeUA = 6) const;
|
||||
|
||||
//@}
|
||||
/// @name Partial Molar Properties of the Solution
|
||||
///
|
||||
/// These properties are handled by the parent class,
|
||||
/// SingleSpeciesTP
|
||||
//@{
|
||||
|
||||
|
||||
//@}
|
||||
/// @name Properties of the Standard State of the Species in the Solution
|
||||
//@{
|
||||
|
||||
//! Get the nondimensional Enthalpy functions for the species
|
||||
//! at their standard states at the current <I>T</I> and <I>P</I> of the solution.
|
||||
/*!
|
||||
* @param hrt Output vector of nondimensional standard state enthalpies.
|
||||
* Length: m_kk.
|
||||
*/
|
||||
virtual void getEnthalpy_RT(doublereal* hrt) const;
|
||||
|
||||
//! Get the array of nondimensional Entropy functions for the
|
||||
//! standard state species at the current <I>T</I> and <I>P</I> of the solution.
|
||||
/*!
|
||||
* @param sr Output vector of nondimensional standard state entropies.
|
||||
* Length: m_kk.
|
||||
*/
|
||||
virtual void getEntropy_R(doublereal* sr) const;
|
||||
|
||||
//! Get the nondimensional Gibbs functions for the species
|
||||
//! in their standard states at the current <I>T</I> and <I>P</I> of the solution.
|
||||
/*!
|
||||
* @param grt Output vector of nondimensional standard state gibbs free energies
|
||||
* Length: m_kk.
|
||||
*/
|
||||
virtual void getGibbs_RT(doublereal* grt) const;
|
||||
|
||||
//! Get the nondimensional Heat Capacities at constant
|
||||
//! pressure for the species standard states
|
||||
//! at the current <I>T</I> and <I>P</I> of the solution
|
||||
/*!
|
||||
* @param cpr Output vector of nondimensional standard state heat capacities
|
||||
* Length: m_kk.
|
||||
*/
|
||||
virtual void getCp_R(doublereal* cpr) const;
|
||||
|
||||
//! Returns the vector of nondimensional Internal Energies of the standard
|
||||
//! state species at the current <I>T</I> and <I>P</I> of the solution
|
||||
/*!
|
||||
* For an incompressible,
|
||||
* stoichiometric substance, the molar internal energy is
|
||||
* independent of pressure. Since the thermodynamic properties
|
||||
* are specified by giving the standard-state enthalpy, the
|
||||
* term \f$ P_{ref} \hat v\f$ is subtracted from the specified reference molar
|
||||
* enthalpy to compute the standard state molar internal energy.
|
||||
*
|
||||
* @param urt output vector of nondimensional standard state
|
||||
* internal energies of the species. Length: m_kk.
|
||||
*/
|
||||
virtual void getIntEnergy_RT(doublereal* urt) const;
|
||||
|
||||
//@}
|
||||
/// @name Thermodynamic Values for the Species Reference States
|
||||
//@{
|
||||
|
||||
//! Returns the vector of nondimensional
|
||||
//! internal Energies of the reference state at the current temperature
|
||||
//! of the solution and the reference pressure for each species.
|
||||
/*!
|
||||
* @param urt Output vector of nondimensional reference state
|
||||
* internal energies of the species.
|
||||
* Length: m_kk
|
||||
*/
|
||||
virtual void getIntEnergy_RT_ref(doublereal *urt) const;
|
||||
|
||||
/*
|
||||
* ---- Critical State Properties
|
||||
*/
|
||||
|
||||
|
||||
/*
|
||||
* ---- Saturation Properties
|
||||
*/
|
||||
|
||||
/*
|
||||
* @internal Initialize. This method is provided to allow
|
||||
* subclasses to perform any initialization required after all
|
||||
* species have been added. For example, it might be used to
|
||||
* resize internal work arrays that must have an entry for
|
||||
* each species. The base class implementation does nothing,
|
||||
* and subclasses that do not require initialization do not
|
||||
* need to overload this method. When importing a CTML phase
|
||||
* description, this method is called just prior to returning
|
||||
* from function importPhase.
|
||||
*
|
||||
* @see importCTML.cpp
|
||||
*/
|
||||
virtual void initThermo();
|
||||
|
||||
|
||||
virtual void initThermoXML(XML_Node& phaseNode, std::string id);
|
||||
|
||||
//! Make the default XML tree
|
||||
/*!
|
||||
* @return Returns a malloced XML tree containing the
|
||||
* default info.
|
||||
*/
|
||||
static XML_Node *makeDefaultXMLTree();
|
||||
|
||||
//! Set the equation of state parameters
|
||||
/*!
|
||||
* @internal
|
||||
* The number and meaning of these depends on the subclass.
|
||||
*
|
||||
* @param n number of parameters
|
||||
* @param c array of \a n coefficients
|
||||
* c[0] = density of phase [ kg/m3 ]
|
||||
*/
|
||||
virtual void setParameters(int n, doublereal * const c);
|
||||
|
||||
//! Get the equation of state parameters in a vector
|
||||
/*!
|
||||
* @internal
|
||||
*
|
||||
* @param n number of parameters
|
||||
* @param c array of \a n coefficients
|
||||
*
|
||||
* For this phase:
|
||||
* - n = 1
|
||||
* - c[0] = density of phase [ kg/m3 ]
|
||||
*/
|
||||
virtual void getParameters(int &n, doublereal * const c) const;
|
||||
|
||||
//! Set equation of state parameter values from XML entries.
|
||||
/*!
|
||||
* This method is called by function importPhase() in
|
||||
* file importCTML.cpp when processing a phase definition in
|
||||
* an input file. It should be overloaded in subclasses to set
|
||||
* any parameters that are specific to that particular phase
|
||||
* model. Note, this method is called before the phase is
|
||||
* initialzed with elements and/or species.
|
||||
*
|
||||
* For this phase, the density of the phase is specified in this block.
|
||||
*
|
||||
* @param eosdata An XML_Node object corresponding to
|
||||
* the "thermo" entry for this phase in the input file.
|
||||
*
|
||||
* eosdata points to the thermo block, and looks like this:
|
||||
*
|
||||
* @verbatim
|
||||
<phase id="stoichsolid" >
|
||||
<thermo model="StoichSubstance">
|
||||
<density units="g/cm3">3.52</density>
|
||||
</thermo>
|
||||
</phase> @endverbatim
|
||||
*
|
||||
*/
|
||||
virtual void setParametersFromXML(const XML_Node& eosdata);
|
||||
|
||||
protected:
|
||||
|
||||
XML_Node *xdef_;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
#endif
|
||||
|
|
@ -543,13 +543,33 @@ namespace Cantera {
|
|||
* @see importCTML.cpp
|
||||
*/
|
||||
void SingleSpeciesTP::initThermo() {
|
||||
|
||||
/*
|
||||
* Check to make sure that there is one and only one species
|
||||
* in this phase.
|
||||
* Make sure there is one and only one species in this phase.
|
||||
*/
|
||||
m_kk = nSpecies();
|
||||
if (m_kk != 1) {
|
||||
err("singleSpeciesTP ERROR m_kk != 1");
|
||||
throw CanteraError("initThermo",
|
||||
"stoichiometric substances may only contain one species.");
|
||||
}
|
||||
doublereal tmin = m_spthermo->minTemp();
|
||||
doublereal tmax = m_spthermo->maxTemp();
|
||||
if (tmin > 0.0) m_tmin = tmin;
|
||||
if (tmax > 0.0) m_tmax = tmax;
|
||||
|
||||
/*
|
||||
* Store the reference pressure in the variables for the class.
|
||||
*/
|
||||
m_p0 = refPressure();
|
||||
|
||||
/*
|
||||
* Resize temporary arrays.
|
||||
*/
|
||||
int leng = 1;
|
||||
m_h0_RT.resize(leng);
|
||||
m_cp0_R.resize(leng);
|
||||
m_s0_R.resize(leng);
|
||||
|
||||
/*
|
||||
* Make sure the species mole fraction is equal to 1.0;
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -59,6 +59,7 @@
|
|||
|
||||
#ifdef WITH_STOICH_SUBSTANCE
|
||||
#include "MineralEQ3.h"
|
||||
#include "MetalSHEelectrons.h"
|
||||
#endif
|
||||
|
||||
//#include "importCTML.h"
|
||||
|
|
@ -93,7 +94,7 @@ namespace Cantera {
|
|||
"PureFluid", "LatticeSolid", "Lattice",
|
||||
"HMW", "IdealSolidSolution", "DebyeHuckel",
|
||||
"IdealMolalSolution", "IdealGasVPSS",
|
||||
"MineralEQ3", "electrodeElectron", "Margules",
|
||||
"MineralEQ3", "MetalSHEelectrons", "Margules",
|
||||
"IonsFromNeutralMolecule"
|
||||
};
|
||||
|
||||
|
|
@ -102,7 +103,7 @@ namespace Cantera {
|
|||
cPureFluid, cLatticeSolid, cLattice,
|
||||
cHMW, cIdealSolidSolnPhase, cDebyeHuckel,
|
||||
cIdealMolalSoln, cVPSS_IdealGas,
|
||||
cMineralEQ3, cElectrodeElectron,
|
||||
cMineralEQ3, cMetalSHEelectrons,
|
||||
cMargulesVPSSTP, cIonsFromNeutral
|
||||
};
|
||||
|
||||
|
|
@ -173,8 +174,8 @@ namespace Cantera {
|
|||
#endif
|
||||
|
||||
#ifdef WITH_STOICH_SUBSTANCE
|
||||
case cElectrodeElectron:
|
||||
th = new electrodeElectron();
|
||||
case cMetalSHEelectrons:
|
||||
th = new MetalSHEelectrons();
|
||||
break;
|
||||
#endif
|
||||
|
||||
|
|
|
|||
|
|
@ -495,6 +495,7 @@ namespace Cantera {
|
|||
|
||||
|
||||
// Parameters for the viscosityWater() function
|
||||
|
||||
//@{
|
||||
const double H[4] = {1.,
|
||||
0.978197,
|
||||
|
|
|
|||
|
|
@ -29,10 +29,10 @@ using std::fabs;
|
|||
* routine, except for internal checks. All calculations here are done
|
||||
* in dimensionless units.
|
||||
*/
|
||||
static const double T_c = 647.096; // Kelvin
|
||||
static const double P_c = 22.064E6; // Pascals
|
||||
static const double Rho_c = 322.; // kg m-3
|
||||
static const double M_water = 18.015268; // kg kmol-1
|
||||
static const doublereal T_c = 647.096; // Kelvin
|
||||
static const doublereal P_c = 22.064E6; // Pascals
|
||||
static const doublereal Rho_c = 322.; // kg m-3
|
||||
static const doublereal M_water = 18.015268; // kg kmol-1
|
||||
|
||||
/*
|
||||
* The added constants were calculated so that u = s = 0
|
||||
|
|
@ -42,7 +42,7 @@ static const double M_water = 18.015268; // kg kmol-1
|
|||
* H didn't turn out to be .611872 J/kg, but .611782 J/kg.
|
||||
* There may be a slight error here somehow.
|
||||
*/
|
||||
static const double ni0[9] = {
|
||||
static const doublereal ni0[9] = {
|
||||
0.0,
|
||||
-8.32044648201 - 0.000000001739715,
|
||||
6.6832105268 + 0.000000000793232,
|
||||
|
|
@ -54,7 +54,7 @@ static const double ni0[9] = {
|
|||
0.24873
|
||||
};
|
||||
|
||||
static const double gammi0[9] = {
|
||||
static const doublereal gammi0[9] = {
|
||||
0.0,
|
||||
0.0,
|
||||
0.0,
|
||||
|
|
@ -241,7 +241,7 @@ static const int tiR[55] = {
|
|||
4 // 54
|
||||
};
|
||||
|
||||
static const double ni[57] = {
|
||||
static const doublereal ni[57] = {
|
||||
+0.0,
|
||||
+0.12533547935523E-1, // 1
|
||||
+0.78957634722828E1, // 2
|
||||
|
|
@ -302,61 +302,61 @@ static const double ni[57] = {
|
|||
};
|
||||
|
||||
|
||||
static const double alphai[3] = {
|
||||
static const doublereal alphai[3] = {
|
||||
+20.,
|
||||
+20.,
|
||||
+20.
|
||||
};
|
||||
|
||||
static const double betai[3] = {
|
||||
static const doublereal betai[3] = {
|
||||
+150.,
|
||||
+150.,
|
||||
+250.
|
||||
};
|
||||
|
||||
static const double gammai[3] = {
|
||||
static const doublereal gammai[3] = {
|
||||
+1.21,
|
||||
+1.21,
|
||||
+1.25
|
||||
};
|
||||
|
||||
static const double epsi[3] = {
|
||||
static const doublereal epsi[3] = {
|
||||
+1.0,
|
||||
+1.0,
|
||||
+1.0
|
||||
};
|
||||
|
||||
static const double ai[2] = {
|
||||
static const doublereal ai[2] = {
|
||||
+3.5,
|
||||
+3.5
|
||||
};
|
||||
|
||||
static const double bi[2] = {
|
||||
static const doublereal bi[2] = {
|
||||
+0.85,
|
||||
+0.95
|
||||
};
|
||||
|
||||
static const double Bi[2] = {
|
||||
static const doublereal Bi[2] = {
|
||||
+0.2,
|
||||
+0.2
|
||||
};
|
||||
|
||||
static const double Ci[2] = {
|
||||
static const doublereal Ci[2] = {
|
||||
+28.0,
|
||||
+32.0
|
||||
};
|
||||
|
||||
static const double Di[2] = {
|
||||
static const doublereal Di[2] = {
|
||||
+700.,
|
||||
+800.
|
||||
};
|
||||
|
||||
static const double Ai[2] = {
|
||||
static const doublereal Ai[2] = {
|
||||
+0.32,
|
||||
+0.32
|
||||
};
|
||||
|
||||
static const double Bbetai[2] = {
|
||||
static const doublereal Bbetai[2] = {
|
||||
+0.3,
|
||||
+0.3
|
||||
};
|
||||
|
|
@ -382,20 +382,20 @@ WaterPropsIAPWSphi::WaterPropsIAPWSphi() :
|
|||
* prints out the result. It's used for conducting the internal
|
||||
* check.
|
||||
*/
|
||||
void WaterPropsIAPWSphi::intCheck(double tau, double delta) {
|
||||
void WaterPropsIAPWSphi::intCheck(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double nau = phi0();
|
||||
double res = phiR();
|
||||
double res_d = phiR_d();
|
||||
double nau_d = phi0_d();
|
||||
double res_dd = phiR_dd();
|
||||
double nau_dd = phi0_dd();
|
||||
double res_t = phiR_t();
|
||||
double nau_t = phi0_t();
|
||||
double res_tt = phiR_tt();
|
||||
double nau_tt = phi0_tt();
|
||||
double res_dt = phiR_dt();
|
||||
double nau_dt = phi0_dt();
|
||||
doublereal nau = phi0();
|
||||
doublereal res = phiR();
|
||||
doublereal res_d = phiR_d();
|
||||
doublereal nau_d = phi0_d();
|
||||
doublereal res_dd = phiR_dd();
|
||||
doublereal nau_dd = phi0_dd();
|
||||
doublereal res_t = phiR_t();
|
||||
doublereal nau_t = phi0_t();
|
||||
doublereal res_tt = phiR_tt();
|
||||
doublereal nau_tt = phi0_tt();
|
||||
doublereal res_dt = phiR_dt();
|
||||
doublereal nau_dt = phi0_dt();
|
||||
|
||||
std::printf("nau = %20.12e\t\tres = %20.12e\n", nau, res);
|
||||
std::printf("nau_d = %20.12e\t\tres_d = %20.12e\n", nau_d, res_d);
|
||||
|
|
@ -406,19 +406,19 @@ void WaterPropsIAPWSphi::intCheck(double tau, double delta) {
|
|||
}
|
||||
|
||||
void WaterPropsIAPWSphi::check1() {
|
||||
double T = 500.;
|
||||
double rho = 838.025;
|
||||
double tau = T_c/T;
|
||||
double delta = rho / Rho_c;
|
||||
doublereal T = 500.;
|
||||
doublereal rho = 838.025;
|
||||
doublereal tau = T_c/T;
|
||||
doublereal delta = rho / Rho_c;
|
||||
printf(" T = 500 K, rho = 838.025 kg m-3\n");
|
||||
intCheck(tau, delta);
|
||||
}
|
||||
|
||||
void WaterPropsIAPWSphi::check2() {
|
||||
double T = 647;
|
||||
double rho = 358.0;
|
||||
double tau = T_c/T;
|
||||
double delta = rho / Rho_c;
|
||||
doublereal T = 647;
|
||||
doublereal rho = 358.0;
|
||||
doublereal tau = T_c/T;
|
||||
doublereal delta = rho / Rho_c;
|
||||
printf(" T = 647 K, rho = 358.0 kg m-3\n");
|
||||
intCheck(tau, delta);
|
||||
}
|
||||
|
|
@ -427,7 +427,7 @@ void WaterPropsIAPWSphi::check2() {
|
|||
* Calculate the polynomials in tau and delta, and store them in static
|
||||
* storage.
|
||||
*/
|
||||
void WaterPropsIAPWSphi::tdpolycalc(double tau, double delta) {
|
||||
void WaterPropsIAPWSphi::tdpolycalc(doublereal tau, doublereal delta) {
|
||||
if ((tau != TAUsave) || 1) {
|
||||
TAUsave = tau;
|
||||
TAUsqrt = sqrt(tau);
|
||||
|
|
@ -449,10 +449,10 @@ void WaterPropsIAPWSphi::tdpolycalc(double tau, double delta) {
|
|||
* Calculate Eqn. 6.5 for phi0, the ideal gas part of the
|
||||
* dimensionless Helmholtz free energy.
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi0() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
double retn = log(delta) + ni0[1] + ni0[2]*tau + ni0[3]*log(tau);
|
||||
doublereal WaterPropsIAPWSphi::phi0() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
doublereal retn = log(delta) + ni0[1] + ni0[2]*tau + ni0[3]*log(tau);
|
||||
|
||||
retn += ni0[4] * log(1.0 - exp(-gammi0[4]*tau));
|
||||
retn += ni0[5] * log(1.0 - exp(-gammi0[5]*tau));
|
||||
|
|
@ -469,16 +469,16 @@ double WaterPropsIAPWSphi::phi0() const {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phiR() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phiR() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
int i, j;
|
||||
|
||||
/*
|
||||
* Write out the first seven polynomials in the expression
|
||||
*/
|
||||
double T375 = pow(tau, 0.375);
|
||||
double val = (ni[1] * delta / TAUsqrt +
|
||||
doublereal T375 = pow(tau, 0.375);
|
||||
doublereal val = (ni[1] * delta / TAUsqrt +
|
||||
ni[2] * delta * TAUsqrt * T375 +
|
||||
ni[3] * delta * tau +
|
||||
ni[4] * DELTAp[2] * TAUsqrt +
|
||||
|
|
@ -497,8 +497,8 @@ double WaterPropsIAPWSphi::phiR() const {
|
|||
*/
|
||||
for (j = 0; j < 3; j++) {
|
||||
i = 52 + j;
|
||||
double dtmp = delta - epsi[j];
|
||||
double ttmp = tau - gammai[j];
|
||||
doublereal dtmp = delta - epsi[j];
|
||||
doublereal ttmp = tau - gammai[j];
|
||||
val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
|
||||
exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
|
||||
}
|
||||
|
|
@ -508,16 +508,16 @@ double WaterPropsIAPWSphi::phiR() const {
|
|||
*/
|
||||
for (j = 0; j < 2; j++) {
|
||||
i = 55 + j;
|
||||
double deltam1 = delta - 1.0;
|
||||
double dtmp2 = deltam1 * deltam1;
|
||||
double atmp = 0.5 / Bbetai[j];
|
||||
double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
double ttmp = tau - 1.0;
|
||||
doublereal deltam1 = delta - 1.0;
|
||||
doublereal dtmp2 = deltam1 * deltam1;
|
||||
doublereal atmp = 0.5 / Bbetai[j];
|
||||
doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
doublereal ttmp = tau - 1.0;
|
||||
|
||||
double triagtmp = pow(triag, bi[j]);
|
||||
doublereal triagtmp = pow(triag, bi[j]);
|
||||
|
||||
double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
val += (ni[i] * triagtmp * delta * phi);
|
||||
}
|
||||
|
||||
|
|
@ -528,11 +528,11 @@ double WaterPropsIAPWSphi::phiR() const {
|
|||
* Calculate the Phi function, which is basically the helmholtz free energy
|
||||
* Eqn. (6.4)
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::phi(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double nau = phi0();
|
||||
double res = phiR();
|
||||
double retn = nau + res;
|
||||
doublereal nau = phi0();
|
||||
doublereal res = phiR();
|
||||
doublereal retn = nau + res;
|
||||
return retn;
|
||||
}
|
||||
|
||||
|
|
@ -544,16 +544,16 @@ double WaterPropsIAPWSphi::phi(double tau, double delta) {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phiR_d() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phiR_d() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
int i, j;
|
||||
|
||||
/*
|
||||
* Write out the first seven polynomials in the expression
|
||||
*/
|
||||
double T375 = pow(tau, 0.375);
|
||||
double val = (ni[1] / TAUsqrt +
|
||||
doublereal T375 = pow(tau, 0.375);
|
||||
doublereal val = (ni[1] / TAUsqrt +
|
||||
ni[2] * TAUsqrt * T375 +
|
||||
ni[3] * tau +
|
||||
ni[4] * 2.0 * delta * TAUsqrt +
|
||||
|
|
@ -573,9 +573,9 @@ double WaterPropsIAPWSphi::phiR_d() const {
|
|||
*/
|
||||
for (j = 0; j < 3; j++) {
|
||||
i = 52 + j;
|
||||
double dtmp = delta - epsi[j];
|
||||
double ttmp = tau - gammai[j];
|
||||
double tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
|
||||
doublereal dtmp = delta - epsi[j];
|
||||
doublereal ttmp = tau - gammai[j];
|
||||
doublereal tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
|
||||
exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
|
||||
val += tmp * (diR[i]/delta - 2.0 * alphai[j] * dtmp);
|
||||
}
|
||||
|
|
@ -585,27 +585,27 @@ double WaterPropsIAPWSphi::phiR_d() const {
|
|||
*/
|
||||
for (j = 0; j < 2; j++) {
|
||||
i = 55 + j;
|
||||
double deltam1 = delta - 1.0;
|
||||
double dtmp2 = deltam1 * deltam1;
|
||||
double atmp = 0.5 / Bbetai[j];
|
||||
double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
double ttmp = tau - 1.0;
|
||||
doublereal deltam1 = delta - 1.0;
|
||||
doublereal dtmp2 = deltam1 * deltam1;
|
||||
doublereal atmp = 0.5 / Bbetai[j];
|
||||
doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
doublereal ttmp = tau - 1.0;
|
||||
|
||||
double triagtmp = pow(triag, bi[j]);
|
||||
double triagtmpm1 = pow(triag, bi[j]-1.0);
|
||||
double atmpM1 = atmp - 1.0;
|
||||
double ptmp = pow(dtmp2,atmpM1);
|
||||
double p2tmp = pow(dtmp2, ai[j]-1.0);
|
||||
double dtriagddelta =
|
||||
doublereal triagtmp = pow(triag, bi[j]);
|
||||
doublereal triagtmpm1 = pow(triag, bi[j]-1.0);
|
||||
doublereal atmpM1 = atmp - 1.0;
|
||||
doublereal ptmp = pow(dtmp2,atmpM1);
|
||||
doublereal p2tmp = pow(dtmp2, ai[j]-1.0);
|
||||
doublereal dtriagddelta =
|
||||
deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
|
||||
2.0*Bi[j]*ai[j]*p2tmp);
|
||||
|
||||
double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
double dphiddelta = -2.0*Ci[j]*deltam1*phi;
|
||||
double dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
|
||||
doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
doublereal dphiddelta = -2.0*Ci[j]*deltam1*phi;
|
||||
doublereal dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
|
||||
|
||||
double tmp = ni[i] * (triagtmp * (phi + delta*dphiddelta) +
|
||||
doublereal tmp = ni[i] * (triagtmp * (phi + delta*dphiddelta) +
|
||||
dtriagtmpddelta * delta * phi);
|
||||
val += tmp;
|
||||
}
|
||||
|
|
@ -620,8 +620,8 @@ double WaterPropsIAPWSphi::phiR_d() const {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi0_d() const {
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phi0_d() const {
|
||||
doublereal delta = DELTAsave;
|
||||
return (1.0/delta);
|
||||
}
|
||||
|
||||
|
|
@ -630,11 +630,11 @@ double WaterPropsIAPWSphi::phi0_d() const {
|
|||
* of helmholtz free energy wrt delta
|
||||
* Eqn. (6.4)
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi_d(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::phi_d(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double nau = phi0_d();
|
||||
double res = phiR_d();
|
||||
double retn = nau + res;
|
||||
doublereal nau = phi0_d();
|
||||
doublereal res = phiR_d();
|
||||
doublereal retn = nau + res;
|
||||
return retn;
|
||||
}
|
||||
|
||||
|
|
@ -645,10 +645,10 @@ double WaterPropsIAPWSphi::phi_d(double tau, double delta) {
|
|||
*
|
||||
* note: this is done so much, we have a seperate routine.
|
||||
*/
|
||||
double WaterPropsIAPWSphi::pressureM_rhoRT(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::pressureM_rhoRT(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double res = phiR_d();
|
||||
double retn = 1.0 + delta * res;
|
||||
doublereal res = phiR_d();
|
||||
doublereal retn = 1.0 + delta * res;
|
||||
return retn;
|
||||
}
|
||||
|
||||
|
|
@ -659,17 +659,17 @@ double WaterPropsIAPWSphi::pressureM_rhoRT(double tau, double delta) {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phiR_dd() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phiR_dd() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
int i, j;
|
||||
double atmp;
|
||||
doublereal atmp;
|
||||
|
||||
/*
|
||||
* Write out the first seven polynomials in the expression
|
||||
*/
|
||||
double T375 = pow(tau, 0.375);
|
||||
double val = (ni[4] * 2.0 * TAUsqrt +
|
||||
doublereal T375 = pow(tau, 0.375);
|
||||
doublereal val = (ni[4] * 2.0 * TAUsqrt +
|
||||
ni[5] * 2.0 * T375 * T375 +
|
||||
ni[6] * 6.0 * delta * T375 +
|
||||
ni[7] * 12.0 * DELTAp[2] * tau);
|
||||
|
|
@ -677,8 +677,8 @@ double WaterPropsIAPWSphi::phiR_dd() const {
|
|||
* Next, do polynomial contributions 8 to 51
|
||||
*/
|
||||
for (i = 8; i <= 51; i++) {
|
||||
double dtmp = DELTAp[ciR[i]];
|
||||
double tmp = ni[i] * exp(-dtmp) * TAUp[tiR[i]];
|
||||
doublereal dtmp = DELTAp[ciR[i]];
|
||||
doublereal tmp = ni[i] * exp(-dtmp) * TAUp[tiR[i]];
|
||||
if (diR[i] == 1) {
|
||||
atmp = 1.0/delta;
|
||||
} else {
|
||||
|
|
@ -694,14 +694,14 @@ double WaterPropsIAPWSphi::phiR_dd() const {
|
|||
*/
|
||||
for (j = 0; j < 3; j++) {
|
||||
i = 52 + j;
|
||||
double dtmp = delta - epsi[j];
|
||||
double ttmp = tau - gammai[j];
|
||||
double tmp = (ni[i] * TAUp[tiR[i]] *
|
||||
doublereal dtmp = delta - epsi[j];
|
||||
doublereal ttmp = tau - gammai[j];
|
||||
doublereal tmp = (ni[i] * TAUp[tiR[i]] *
|
||||
exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
|
||||
double deltmp = DELTAp[diR[i]];
|
||||
double deltmpM1 = deltmp/delta;
|
||||
double deltmpM2 = deltmpM1 / delta;
|
||||
double d2tmp = dtmp * dtmp;
|
||||
doublereal deltmp = DELTAp[diR[i]];
|
||||
doublereal deltmpM1 = deltmp/delta;
|
||||
doublereal deltmpM2 = deltmpM1 / delta;
|
||||
doublereal d2tmp = dtmp * dtmp;
|
||||
|
||||
val += tmp * (-2.0*alphai[j]*deltmp +
|
||||
4.0 * alphai[j] * alphai[j] * deltmp * d2tmp -
|
||||
|
|
@ -714,41 +714,41 @@ double WaterPropsIAPWSphi::phiR_dd() const {
|
|||
*/
|
||||
for (j = 0; j < 2; j++) {
|
||||
i = 55 + j;
|
||||
double deltam1 = delta - 1.0;
|
||||
double dtmp2 = deltam1 * deltam1;
|
||||
doublereal deltam1 = delta - 1.0;
|
||||
doublereal dtmp2 = deltam1 * deltam1;
|
||||
atmp = 0.5 / Bbetai[j];
|
||||
double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
double ttmp = tau - 1.0;
|
||||
doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
doublereal ttmp = tau - 1.0;
|
||||
|
||||
double triagtmp = pow(triag, bi[j]);
|
||||
double triagtmpm1 = pow(triag, bi[j]-1.0);
|
||||
double atmpM1 = atmp - 1.0;
|
||||
double ptmp = pow(dtmp2,atmpM1);
|
||||
double p2tmp = pow(dtmp2, ai[j]-1.0);
|
||||
double dtriagddelta =
|
||||
doublereal triagtmp = pow(triag, bi[j]);
|
||||
doublereal triagtmpm1 = pow(triag, bi[j]-1.0);
|
||||
doublereal atmpM1 = atmp - 1.0;
|
||||
doublereal ptmp = pow(dtmp2,atmpM1);
|
||||
doublereal p2tmp = pow(dtmp2, ai[j]-1.0);
|
||||
doublereal dtriagddelta =
|
||||
deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
|
||||
2.0*Bi[j]*ai[j]*p2tmp);
|
||||
|
||||
double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
double dphiddelta = -2.0*Ci[j]*deltam1*phi;
|
||||
double dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
|
||||
doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
doublereal dphiddelta = -2.0*Ci[j]*deltam1*phi;
|
||||
doublereal dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
|
||||
|
||||
|
||||
double d2phiddelta2 = 2.0 * Ci[j] * phi * (2.0*Ci[j]*dtmp2 - 1.0);
|
||||
doublereal d2phiddelta2 = 2.0 * Ci[j] * phi * (2.0*Ci[j]*dtmp2 - 1.0);
|
||||
|
||||
double pptmp = ptmp / dtmp2;
|
||||
double d2triagddelta2 = dtriagddelta / deltam1;
|
||||
doublereal pptmp = ptmp / dtmp2;
|
||||
doublereal d2triagddelta2 = dtriagddelta / deltam1;
|
||||
d2triagddelta2 +=
|
||||
dtmp2 *(4.0*Bi[j]*ai[j]*(ai[j]-1.0)*pow(dtmp2,ai[j]-2.0) +
|
||||
2.0*Ai[j]*Ai[j]/(Bbetai[j]*Bbetai[j])*ptmp*ptmp +
|
||||
Ai[j]*theta*4.0/Bbetai[j]*(atmp-1.0)*pptmp);
|
||||
|
||||
double d2triagtmpd2delta =
|
||||
doublereal d2triagtmpd2delta =
|
||||
bi[j] * (triagtmpm1 * d2triagddelta2 +
|
||||
(bi[j]-1.0)*triagtmpm1/triag*dtriagddelta*dtriagddelta);
|
||||
|
||||
double ctmp = (triagtmp * (2.0*dphiddelta + delta*d2phiddelta2) +
|
||||
doublereal ctmp = (triagtmp * (2.0*dphiddelta + delta*d2phiddelta2) +
|
||||
2.0*dtriagtmpddelta*(phi + delta * dphiddelta) +
|
||||
d2triagtmpd2delta * delta * phi);
|
||||
|
||||
|
|
@ -765,8 +765,8 @@ double WaterPropsIAPWSphi::phiR_dd() const {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi0_dd() const {
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phi0_dd() const {
|
||||
doublereal delta = DELTAsave;
|
||||
return (-1.0/(delta*delta));
|
||||
}
|
||||
|
||||
|
|
@ -775,36 +775,36 @@ double WaterPropsIAPWSphi::phi0_dd() const {
|
|||
* of helmholtz free energy wrt delta
|
||||
* Eqn. (6.4)
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi_dd(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::phi_dd(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double nau = phi0_dd();
|
||||
double res = phiR_dd();
|
||||
double retn = nau + res;
|
||||
doublereal nau = phi0_dd();
|
||||
doublereal res = phiR_dd();
|
||||
doublereal retn = nau + res;
|
||||
return retn;
|
||||
}
|
||||
|
||||
double WaterPropsIAPWSphi::dimdpdrho(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::dimdpdrho(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double res1 = phiR_d();
|
||||
double res2 = phiR_dd();
|
||||
double retn = 1.0 + delta * (2.0*res1 + delta*res2);
|
||||
doublereal res1 = phiR_d();
|
||||
doublereal res2 = phiR_dd();
|
||||
doublereal retn = 1.0 + delta * (2.0*res1 + delta*res2);
|
||||
return retn;
|
||||
}
|
||||
|
||||
double WaterPropsIAPWSphi::dimdpdT(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::dimdpdT(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double res1 = phiR_d();
|
||||
double res2 = phiR_dt();
|
||||
double retn = (1.0 + delta * res1) - tau * delta * (res2);
|
||||
doublereal res1 = phiR_d();
|
||||
doublereal res2 = phiR_dt();
|
||||
doublereal retn = (1.0 + delta * res1) - tau * delta * (res2);
|
||||
return retn;
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate d_phi0/d(tau)
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi0_t() const {
|
||||
double tau = TAUsave;
|
||||
double retn = ni0[2] + ni0[3]/tau;;
|
||||
doublereal WaterPropsIAPWSphi::phi0_t() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal retn = ni0[2] + ni0[3]/tau;;
|
||||
retn += (ni0[4] * gammi0[4] * (1.0/(1.0 - exp(-gammi0[4]*tau)) - 1.0));
|
||||
retn += (ni0[5] * gammi0[5] * (1.0/(1.0 - exp(-gammi0[5]*tau)) - 1.0));
|
||||
retn += (ni0[6] * gammi0[6] * (1.0/(1.0 - exp(-gammi0[6]*tau)) - 1.0));
|
||||
|
|
@ -820,17 +820,17 @@ double WaterPropsIAPWSphi::phi0_t() const {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phiR_t() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phiR_t() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
int i, j;
|
||||
double atmp, tmp;
|
||||
doublereal atmp, tmp;
|
||||
|
||||
/*
|
||||
* Write out the first seven polynomials in the expression
|
||||
*/
|
||||
double T375 = pow(tau, 0.375);
|
||||
double val = ((-0.5) *ni[1] * delta / TAUsqrt / tau +
|
||||
doublereal T375 = pow(tau, 0.375);
|
||||
doublereal val = ((-0.5) *ni[1] * delta / TAUsqrt / tau +
|
||||
ni[2] * delta * 0.875 / TAUsqrt * T375 +
|
||||
ni[3] * delta +
|
||||
ni[4] * DELTAp[2] * 0.5 / TAUsqrt +
|
||||
|
|
@ -850,8 +850,8 @@ double WaterPropsIAPWSphi::phiR_t() const {
|
|||
*/
|
||||
for (j = 0; j < 3; j++) {
|
||||
i = 52 + j;
|
||||
double dtmp = delta - epsi[j];
|
||||
double ttmp = tau - gammai[j];
|
||||
doublereal dtmp = delta - epsi[j];
|
||||
doublereal ttmp = tau - gammai[j];
|
||||
tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
|
||||
exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
|
||||
val += tmp *(tiR[i]/tau - 2.0 * betai[j]*ttmp);
|
||||
|
|
@ -862,21 +862,21 @@ double WaterPropsIAPWSphi::phiR_t() const {
|
|||
*/
|
||||
for (j = 0; j < 2; j++) {
|
||||
i = 55 + j;
|
||||
double deltam1 = delta - 1.0;
|
||||
double dtmp2 = deltam1 * deltam1;
|
||||
doublereal deltam1 = delta - 1.0;
|
||||
doublereal dtmp2 = deltam1 * deltam1;
|
||||
atmp = 0.5 / Bbetai[j];
|
||||
double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
double ttmp = tau - 1.0;
|
||||
doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
doublereal ttmp = tau - 1.0;
|
||||
|
||||
double triagtmp = pow(triag, bi[j]);
|
||||
doublereal triagtmp = pow(triag, bi[j]);
|
||||
|
||||
double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
|
||||
|
||||
double dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
|
||||
doublereal dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
|
||||
|
||||
double dphidtau = - 2.0 * Di[j] * ttmp * phi;
|
||||
doublereal dphidtau = - 2.0 * Di[j] * ttmp * phi;
|
||||
|
||||
val += ni[i] * delta * (dtriagtmpdtau * phi + triagtmp * dphidtau);
|
||||
}
|
||||
|
|
@ -889,21 +889,21 @@ double WaterPropsIAPWSphi::phiR_t() const {
|
|||
* of helmholtz free energy wrt tau
|
||||
* Eqn. (6.4)
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi_t(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::phi_t(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double nau = phi0_t();
|
||||
double res = phiR_t();
|
||||
double retn = nau + res;
|
||||
doublereal nau = phi0_t();
|
||||
doublereal res = phiR_t();
|
||||
doublereal retn = nau + res;
|
||||
return retn;
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate d2_phi0/dtau2
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi0_tt() const {
|
||||
double tau = TAUsave;
|
||||
double tmp, itmp;
|
||||
double retn = - ni0[3]/(tau * tau);
|
||||
doublereal WaterPropsIAPWSphi::phi0_tt() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal tmp, itmp;
|
||||
doublereal retn = - ni0[3]/(tau * tau);
|
||||
for (int i = 4; i <= 8; i++) {
|
||||
tmp = exp(-gammi0[i]*tau);
|
||||
itmp = 1.0 - tmp;
|
||||
|
|
@ -919,17 +919,17 @@ double WaterPropsIAPWSphi::phi0_tt() const {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phiR_tt() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phiR_tt() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
int i, j;
|
||||
double atmp, tmp;
|
||||
doublereal atmp, tmp;
|
||||
|
||||
/*
|
||||
* Write out the first seven polynomials in the expression
|
||||
*/
|
||||
double T375 = pow(tau, 0.375);
|
||||
double val = ((-0.5) * (-1.5) * ni[1] * delta / (TAUsqrt * tau * tau) +
|
||||
doublereal T375 = pow(tau, 0.375);
|
||||
doublereal val = ((-0.5) * (-1.5) * ni[1] * delta / (TAUsqrt * tau * tau) +
|
||||
ni[2] * delta * 0.875 * (-0.125) * T375 / (TAUsqrt * tau) +
|
||||
ni[4] * DELTAp[2] * 0.5 * (-0.5)/ (TAUsqrt * tau) +
|
||||
ni[5] * DELTAp[2] * 0.75 *(-0.25) * T375 * T375 / (tau * tau) +
|
||||
|
|
@ -949,8 +949,8 @@ double WaterPropsIAPWSphi::phiR_tt() const {
|
|||
*/
|
||||
for (j = 0; j < 3; j++) {
|
||||
i = 52 + j;
|
||||
double dtmp = delta - epsi[j];
|
||||
double ttmp = tau - gammai[j];
|
||||
doublereal dtmp = delta - epsi[j];
|
||||
doublereal ttmp = tau - gammai[j];
|
||||
tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
|
||||
exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
|
||||
atmp = tiR[i]/tau - 2.0 * betai[j]*ttmp;
|
||||
|
|
@ -962,28 +962,28 @@ double WaterPropsIAPWSphi::phiR_tt() const {
|
|||
*/
|
||||
for (j = 0; j < 2; j++) {
|
||||
i = 55 + j;
|
||||
double deltam1 = delta - 1.0;
|
||||
double dtmp2 = deltam1 * deltam1;
|
||||
doublereal deltam1 = delta - 1.0;
|
||||
doublereal dtmp2 = deltam1 * deltam1;
|
||||
atmp = 0.5 / Bbetai[j];
|
||||
double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
double ttmp = tau - 1.0;
|
||||
doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
doublereal ttmp = tau - 1.0;
|
||||
|
||||
double triagtmp = pow(triag, bi[j]);
|
||||
double triagtmpM1 = triagtmp / triag;
|
||||
doublereal triagtmp = pow(triag, bi[j]);
|
||||
doublereal triagtmpM1 = triagtmp / triag;
|
||||
|
||||
double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
|
||||
|
||||
double dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
|
||||
doublereal dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
|
||||
|
||||
double dphidtau = - 2.0 * Di[j] * ttmp * phi;
|
||||
doublereal dphidtau = - 2.0 * Di[j] * ttmp * phi;
|
||||
|
||||
double d2triagtmpdtau2 =
|
||||
doublereal d2triagtmpdtau2 =
|
||||
(2 * bi[j] * triagtmpM1 +
|
||||
4 * theta * theta * bi[j] * (bi[j]-1.0) * triagtmpM1 / triag);
|
||||
|
||||
double d2phidtau2 = 2.0*Di[j]*phi *(2.0*Di[j]*ttmp*ttmp - 1.0);
|
||||
doublereal d2phidtau2 = 2.0*Di[j]*phi *(2.0*Di[j]*ttmp*ttmp - 1.0);
|
||||
|
||||
tmp = (d2triagtmpdtau2 * phi +
|
||||
2 * dtriagtmpdtau * dphidtau +
|
||||
|
|
@ -999,18 +999,18 @@ double WaterPropsIAPWSphi::phiR_tt() const {
|
|||
* of helmholtz free energy wrt tau
|
||||
* Eqn. (6.4)
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi_tt(double tau, double delta) {
|
||||
doublereal WaterPropsIAPWSphi::phi_tt(doublereal tau, doublereal delta) {
|
||||
tdpolycalc(tau, delta);
|
||||
double nau = phi0_tt();
|
||||
double res = phiR_tt();
|
||||
double retn = nau + res;
|
||||
doublereal nau = phi0_tt();
|
||||
doublereal res = phiR_tt();
|
||||
doublereal retn = nau + res;
|
||||
return retn;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate d2_phi0/dtauddelta
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phi0_dt() const {
|
||||
doublereal WaterPropsIAPWSphi::phi0_dt() const {
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
|
|
@ -1021,16 +1021,16 @@ double WaterPropsIAPWSphi::phi0_dt() const {
|
|||
* tau = dimensionless temperature
|
||||
* delta = dimensionless pressure
|
||||
*/
|
||||
double WaterPropsIAPWSphi::phiR_dt() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
doublereal WaterPropsIAPWSphi::phiR_dt() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
int i, j;
|
||||
double tmp;
|
||||
doublereal tmp;
|
||||
/*
|
||||
* Write out the first seven polynomials in the expression
|
||||
*/
|
||||
double T375 = pow(tau, 0.375);
|
||||
double val = (ni[1] * (-0.5) / (TAUsqrt * tau) +
|
||||
doublereal T375 = pow(tau, 0.375);
|
||||
doublereal val = (ni[1] * (-0.5) / (TAUsqrt * tau) +
|
||||
ni[2] * (0.875) * T375 / TAUsqrt +
|
||||
ni[3] +
|
||||
ni[4] * 2.0 * delta * (0.5) / TAUsqrt +
|
||||
|
|
@ -1051,8 +1051,8 @@ double WaterPropsIAPWSphi::phiR_dt() const {
|
|||
*/
|
||||
for (j = 0; j < 3; j++) {
|
||||
i = 52 + j;
|
||||
double dtmp = delta - epsi[j];
|
||||
double ttmp = tau - gammai[j];
|
||||
doublereal dtmp = delta - epsi[j];
|
||||
doublereal ttmp = tau - gammai[j];
|
||||
tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] *
|
||||
exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
|
||||
val += tmp * ((diR[i]/delta - 2.0 * alphai[j] * dtmp) *
|
||||
|
|
@ -1064,39 +1064,39 @@ double WaterPropsIAPWSphi::phiR_dt() const {
|
|||
*/
|
||||
for (j = 0; j < 2; j++) {
|
||||
i = 55 + j;
|
||||
double deltam1 = delta - 1.0;
|
||||
double dtmp2 = deltam1 * deltam1;
|
||||
double atmp = 0.5 / Bbetai[j];
|
||||
double theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
double triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
double ttmp = tau - 1.0;
|
||||
doublereal deltam1 = delta - 1.0;
|
||||
doublereal dtmp2 = deltam1 * deltam1;
|
||||
doublereal atmp = 0.5 / Bbetai[j];
|
||||
doublereal theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
|
||||
doublereal triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
|
||||
doublereal ttmp = tau - 1.0;
|
||||
|
||||
double triagtmp = pow(triag, bi[j]);
|
||||
double triagtmpm1 = pow(triag, bi[j]-1.0);
|
||||
double atmpM1 = atmp - 1.0;
|
||||
double ptmp = pow(dtmp2,atmpM1);
|
||||
double p2tmp = pow(dtmp2, ai[j]-1.0);
|
||||
double dtriagddelta =
|
||||
doublereal triagtmp = pow(triag, bi[j]);
|
||||
doublereal triagtmpm1 = pow(triag, bi[j]-1.0);
|
||||
doublereal atmpM1 = atmp - 1.0;
|
||||
doublereal ptmp = pow(dtmp2,atmpM1);
|
||||
doublereal p2tmp = pow(dtmp2, ai[j]-1.0);
|
||||
doublereal dtriagddelta =
|
||||
deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
|
||||
2.0*Bi[j]*ai[j]*p2tmp);
|
||||
|
||||
double phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
double dphiddelta = -2.0*Ci[j]*deltam1*phi;
|
||||
double dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
|
||||
doublereal phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
|
||||
doublereal dphiddelta = -2.0*Ci[j]*deltam1*phi;
|
||||
doublereal dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;
|
||||
|
||||
|
||||
double dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
|
||||
doublereal dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;
|
||||
|
||||
double dphidtau = - 2.0 * Di[j] * ttmp * phi;
|
||||
doublereal dphidtau = - 2.0 * Di[j] * ttmp * phi;
|
||||
|
||||
double d2phiddeltadtau = 4.0 * Ci[j] * Di[j] * deltam1 * ttmp * phi;
|
||||
doublereal d2phiddeltadtau = 4.0 * Ci[j] * Di[j] * deltam1 * ttmp * phi;
|
||||
|
||||
double d2triagtmpddeltadtau =
|
||||
doublereal d2triagtmpddeltadtau =
|
||||
( -Ai[j] * bi[j] * 2.0 / Bbetai[j] * triagtmpm1 * deltam1 * ptmp
|
||||
-2.0 * theta * bi[j] * (bi[j] - 1.0) * triagtmpm1 / triag * dtriagddelta);
|
||||
|
||||
|
||||
double tmp = ni[i] * (triagtmp * (dphidtau + delta*d2phiddeltadtau) +
|
||||
doublereal tmp = ni[i] * (triagtmp * (dphidtau + delta*d2phiddeltadtau) +
|
||||
delta * dtriagtmpddelta * dphidtau +
|
||||
dtriagtmpdtau * (phi + delta * dphiddelta) +
|
||||
d2triagtmpddeltadtau * delta * phi);
|
||||
|
|
@ -1113,30 +1113,30 @@ double WaterPropsIAPWSphi::phiR_dt() const {
|
|||
* critical point.
|
||||
*
|
||||
*/
|
||||
double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) {
|
||||
double dd = deltaGuess;
|
||||
doublereal WaterPropsIAPWSphi::dfind(doublereal p_red, doublereal tau, doublereal deltaGuess) {
|
||||
doublereal dd = deltaGuess;
|
||||
bool conv = false;
|
||||
double deldd = dd;
|
||||
double pcheck = 1.0E-30 + 1.0E-8 * p_red;
|
||||
doublereal deldd = dd;
|
||||
doublereal pcheck = 1.0E-30 + 1.0E-8 * p_red;
|
||||
for (int n = 0; n < 200; n++) {
|
||||
/*
|
||||
* Calculate the internal polynomials, and then calculate the
|
||||
* phi deriv functions needed by this routine.
|
||||
*/
|
||||
tdpolycalc(tau, dd);
|
||||
double q1 = phiR_d();
|
||||
double q2 = phiR_dd();
|
||||
doublereal q1 = phiR_d();
|
||||
doublereal q2 = phiR_dd();
|
||||
|
||||
/*
|
||||
* Calculate the predicted reduced pressure, pred0, based on the
|
||||
* current tau and dd.
|
||||
*/
|
||||
double pred0 = dd + dd * dd * q1;
|
||||
doublereal pred0 = dd + dd * dd * q1;
|
||||
/*
|
||||
* Calculate the derivative of the predicted reduced pressure
|
||||
* wrt the reduced density, dd, This is dpddelta
|
||||
*/
|
||||
double dpddelta = 1.0 + 2.0 * dd * q1 + dd * dd * q2;
|
||||
doublereal dpddelta = 1.0 + 2.0 * dd * q1 + dd * dd * q2;
|
||||
/*
|
||||
* If dpddelta is negative, then we are in the middle of the
|
||||
* 2 phase region, beyond the stability curve. We need to adjust
|
||||
|
|
@ -1158,7 +1158,7 @@ double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) {
|
|||
/*
|
||||
* Dampen and crop the update
|
||||
*/
|
||||
double dpdx = dpddelta;
|
||||
doublereal dpdx = dpddelta;
|
||||
if (n < 10) {
|
||||
dpdx = dpddelta * 1.1;
|
||||
}
|
||||
|
|
@ -1200,74 +1200,74 @@ double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) {
|
|||
/**
|
||||
* Calculate the dimensionless gibbs free energy g/RT.
|
||||
*/
|
||||
double WaterPropsIAPWSphi::gibbs_RT() const {
|
||||
double delta = DELTAsave;
|
||||
double rd = phiR_d();
|
||||
double g = 1.0 + phi0() + phiR() + delta * rd;
|
||||
doublereal WaterPropsIAPWSphi::gibbs_RT() const {
|
||||
doublereal delta = DELTAsave;
|
||||
doublereal rd = phiR_d();
|
||||
doublereal g = 1.0 + phi0() + phiR() + delta * rd;
|
||||
return g;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate the dimensionless enthalpy h/RT.
|
||||
*/
|
||||
double WaterPropsIAPWSphi::enthalpy_RT() const {
|
||||
double delta = DELTAsave;
|
||||
double tau = TAUsave;
|
||||
double rd = phiR_d();
|
||||
double nt = phi0_t();
|
||||
double rt = phiR_t();
|
||||
double hRT = 1.0 + tau * (nt + rt) + delta * rd;
|
||||
doublereal WaterPropsIAPWSphi::enthalpy_RT() const {
|
||||
doublereal delta = DELTAsave;
|
||||
doublereal tau = TAUsave;
|
||||
doublereal rd = phiR_d();
|
||||
doublereal nt = phi0_t();
|
||||
doublereal rt = phiR_t();
|
||||
doublereal hRT = 1.0 + tau * (nt + rt) + delta * rd;
|
||||
return hRT;
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the dimensionless entropy s/R.
|
||||
*/
|
||||
double WaterPropsIAPWSphi::entropy_R() const {
|
||||
double tau = TAUsave;
|
||||
double nt = phi0_t();
|
||||
double rt = phiR_t();
|
||||
double p0 = phi0();
|
||||
double pR = phiR();
|
||||
double sR = tau * (nt + rt) - p0 - pR;
|
||||
doublereal WaterPropsIAPWSphi::entropy_R() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal nt = phi0_t();
|
||||
doublereal rt = phiR_t();
|
||||
doublereal p0 = phi0();
|
||||
doublereal pR = phiR();
|
||||
doublereal sR = tau * (nt + rt) - p0 - pR;
|
||||
return sR;
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the dimensionless internal energy, u/RT.
|
||||
*/
|
||||
double WaterPropsIAPWSphi::intEnergy_RT() const {
|
||||
double tau = TAUsave;
|
||||
double nt = phi0_t();
|
||||
double rt = phiR_t();
|
||||
double uR = tau * (nt + rt);
|
||||
doublereal WaterPropsIAPWSphi::intEnergy_RT() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal nt = phi0_t();
|
||||
doublereal rt = phiR_t();
|
||||
doublereal uR = tau * (nt + rt);
|
||||
return uR;
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the dimensionless constant volume Heat Capacity, Cv/R
|
||||
*/
|
||||
double WaterPropsIAPWSphi::cv_R() const {
|
||||
double tau = TAUsave;
|
||||
double ntt = phi0_tt();
|
||||
double rtt = phiR_tt();
|
||||
double cvR = - tau * tau * (ntt + rtt);
|
||||
doublereal WaterPropsIAPWSphi::cv_R() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal ntt = phi0_tt();
|
||||
doublereal rtt = phiR_tt();
|
||||
doublereal cvR = - tau * tau * (ntt + rtt);
|
||||
return cvR;
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the dimensionless constant pressure Heat Capacity, Cp/R
|
||||
*/
|
||||
double WaterPropsIAPWSphi::cp_R() const {
|
||||
double tau = TAUsave;
|
||||
double delta = DELTAsave;
|
||||
double cvR = cv_R();
|
||||
//double nd = phi0_d();
|
||||
double rd = phiR_d();
|
||||
double rdd = phiR_dd();
|
||||
double rdt = phiR_dt();
|
||||
double num = (1.0 + delta * rd - delta * tau * rdt);
|
||||
double cpR = cvR + (num * num /
|
||||
doublereal WaterPropsIAPWSphi::cp_R() const {
|
||||
doublereal tau = TAUsave;
|
||||
doublereal delta = DELTAsave;
|
||||
doublereal cvR = cv_R();
|
||||
//doublereal nd = phi0_d();
|
||||
doublereal rd = phiR_d();
|
||||
doublereal rdd = phiR_dd();
|
||||
doublereal rdt = phiR_dt();
|
||||
doublereal num = (1.0 + delta * rd - delta * tau * rdt);
|
||||
doublereal cpR = cvR + (num * num /
|
||||
(1.0 + 2.0 * delta * rd + delta * delta * rdd));
|
||||
return cpR;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -45,35 +45,35 @@ public:
|
|||
* @param tau Dimensionless temperature = T_c/T
|
||||
* @param delta Dimensionless density = delta = rho / Rho_c
|
||||
*/
|
||||
double phi(double tau, double delta);
|
||||
doublereal phi(doublereal tau, doublereal 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);
|
||||
doublereal phi_d(doublereal tau, doublereal 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);
|
||||
double phi_dd(doublereal tau, doublereal 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);
|
||||
doublereal phi_t(doublereal tau, doublereal 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);
|
||||
doublereal phi_tt(doublereal tau, doublereal delta);
|
||||
|
||||
//! Internal check # 1
|
||||
void check1();
|
||||
|
|
@ -91,7 +91,7 @@ public:
|
|||
*
|
||||
* note: this is done so much, we have a seperate routine.
|
||||
*/
|
||||
double pressureM_rhoRT(double tau, double delta);
|
||||
doublereal pressureM_rhoRT(doublereal tau, doublereal delta);
|
||||
|
||||
//! Dimensionless derivative of p wrt rho at constant T
|
||||
/*!
|
||||
|
|
@ -101,7 +101,7 @@ public:
|
|||
* @param tau Dimensionless temperature = T_c/T
|
||||
* @param delta Dimensionless density = delta = rho / Rho_c
|
||||
*/
|
||||
double dimdpdrho(double tau, double delta);
|
||||
doublereal dimdpdrho(doublereal tau, doublereal delta);
|
||||
|
||||
//! Dimensionless derivative of p wrt T at constant rho
|
||||
/*!
|
||||
|
|
@ -111,7 +111,7 @@ public:
|
|||
* @param tau Dimensionless temperature = T_c/T
|
||||
* @param delta Dimensionless density = delta = rho / Rho_c
|
||||
*/
|
||||
double dimdpdT(double tau, double delta);
|
||||
doublereal dimdpdT(doublereal tau, doublereal delta);
|
||||
|
||||
/**
|
||||
* This program computes the reduced density, given the reduced pressure
|
||||
|
|
@ -126,37 +126,37 @@ public:
|
|||
* @return
|
||||
* Returns the dimensionless density.
|
||||
*/
|
||||
double dfind(double p_red, double tau, double deltaGuess);
|
||||
doublereal dfind(doublereal p_red, doublereal tau, doublereal deltaGuess);
|
||||
|
||||
/**
|
||||
* Calculate the dimensionless gibbs free energy
|
||||
*/
|
||||
double gibbs_RT() const;
|
||||
doublereal gibbs_RT() const;
|
||||
|
||||
/**
|
||||
* Calculate the dimensionless enthalpy, h/RT
|
||||
*/
|
||||
double enthalpy_RT() const;
|
||||
doublereal enthalpy_RT() const;
|
||||
|
||||
/**
|
||||
* Calculate the dimensionless entropy, s/R
|
||||
*/
|
||||
double entropy_R() const;
|
||||
doublereal entropy_R() const;
|
||||
|
||||
/**
|
||||
* Calculate the dimensionless internal energy, u/RT
|
||||
*/
|
||||
double intEnergy_RT() const;
|
||||
doublereal intEnergy_RT() const;
|
||||
|
||||
/**
|
||||
* Calculate the dimensionless constant volume heat capacity, Cv/R
|
||||
*/
|
||||
double cv_R() const;
|
||||
doublereal cv_R() const;
|
||||
|
||||
/**
|
||||
* Calculate the dimensionless constant pressure heat capacity, Cv/R
|
||||
*/
|
||||
double cp_R() const;
|
||||
doublereal cp_R() const;
|
||||
|
||||
|
||||
//! Calculates internal polynomials in tau and delta.
|
||||
|
|
@ -167,35 +167,35 @@ public:
|
|||
* @param tau Dimensionless temperature = T_c/T
|
||||
* @param delta Dimensionless density = delta = rho / Rho_c
|
||||
*/
|
||||
void tdpolycalc(double tau, double delta);
|
||||
void tdpolycalc(doublereal tau, doublereal delta);
|
||||
|
||||
//! Return the value of phiR(), res
|
||||
double phiR() const;
|
||||
doublereal phiR() const;
|
||||
|
||||
private:
|
||||
|
||||
//! nau calculation
|
||||
double phi0() const;
|
||||
doublereal phi0() const;
|
||||
//! calculation of d_phiR/d_d
|
||||
double phiR_d() const;
|
||||
doublereal phiR_d() const;
|
||||
//! calculation of d_nau/d_d
|
||||
double phi0_d() const;
|
||||
doublereal phi0_d() const;
|
||||
//! calculation of d2_res/d_dd
|
||||
double phiR_dd() const;
|
||||
doublereal phiR_dd() const;
|
||||
//! calculation of d2_nau/d_dd
|
||||
double phi0_dd() const;
|
||||
doublereal phi0_dd() const;
|
||||
//! calculation of d_nau/d_t
|
||||
double phi0_t() const;
|
||||
doublereal phi0_t() const;
|
||||
//! calculation of d_res/d_t
|
||||
double phiR_t() const;
|
||||
doublereal phiR_t() const;
|
||||
//! calculation of d2_res/d_tt
|
||||
double phiR_tt() const;
|
||||
doublereal phiR_tt() const;
|
||||
//! calculation of d2_nau/d_tt
|
||||
double phi0_tt() const;
|
||||
doublereal phi0_tt() const;
|
||||
//! calculation of d2_res/d_dt
|
||||
double phiR_dt() const;
|
||||
doublereal phiR_dt() const;
|
||||
//! calculation of d2_nau/d_dt
|
||||
double phi0_dt() const;
|
||||
doublereal phi0_dt() const;
|
||||
|
||||
/**
|
||||
* intCheck() calculates all of the functions at a one point and
|
||||
|
|
@ -205,23 +205,23 @@ private:
|
|||
* @param tau Dimensionless temperature = T_c/T
|
||||
* @param delta Dimensionless density = delta = rho / Rho_c
|
||||
*/
|
||||
void intCheck(double tau, double delta);
|
||||
void intCheck(doublereal tau, doublereal delta);
|
||||
|
||||
private:
|
||||
|
||||
//! Value of internally calculated polynomials of powers of TAU
|
||||
double TAUp[52];
|
||||
doublereal TAUp[52];
|
||||
|
||||
//! Value of internally calculated polynomials of powers of delta
|
||||
double DELTAp[16];
|
||||
doublereal DELTAp[16];
|
||||
|
||||
//! Last tau that was used to calculate polynomials
|
||||
double TAUsave;
|
||||
doublereal TAUsave;
|
||||
|
||||
//! sqrt of TAU
|
||||
double TAUsqrt;
|
||||
doublereal TAUsqrt;
|
||||
|
||||
//! Last delta that was used to calculate polynomials
|
||||
double DELTAsave;
|
||||
doublereal DELTAsave;
|
||||
};
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -44,7 +44,7 @@ namespace Cantera {
|
|||
const int cSemiconductor = 7;
|
||||
|
||||
const int cMineralEQ3 = 8; // MineralEQ3 in MineralEQ3.h
|
||||
const int cElectrodeElectron = 9; // electrodeElectron
|
||||
const int cMetalSHEelectrons = 9; // SHE electrode electrons
|
||||
|
||||
const int cLatticeSolid = 20; // LatticeSolidPhase.h
|
||||
const int cLattice = 21;
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue