From 6b1a1ac27f4fe9c9df5f09fbccbf3d1efbbec4fa Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 11 Dec 2009 15:27:02 +0000 Subject: [PATCH] Merged changes from the trunk. --- Cantera/src/thermo/IdealGasPhase.h | 2 +- Cantera/src/thermo/Makefile.in | 4 +- Cantera/src/thermo/MetalSHEelectrons.cpp | 546 +++++++++++++++++++++ Cantera/src/thermo/MetalSHEelectrons.h | 565 ++++++++++++++++++++++ Cantera/src/thermo/SingleSpeciesTP.cpp | 26 +- Cantera/src/thermo/ThermoFactory.cpp | 9 +- Cantera/src/thermo/WaterProps.cpp | 1 + Cantera/src/thermo/WaterPropsIAPWSphi.cpp | 544 ++++++++++----------- Cantera/src/thermo/WaterPropsIAPWSphi.h | 68 +-- Cantera/src/thermo/mix_defs.h | 2 +- 10 files changed, 1450 insertions(+), 317 deletions(-) create mode 100644 Cantera/src/thermo/MetalSHEelectrons.cpp create mode 100644 Cantera/src/thermo/MetalSHEelectrons.h diff --git a/Cantera/src/thermo/IdealGasPhase.h b/Cantera/src/thermo/IdealGasPhase.h index cc676a8a4..b83cec970 100644 --- a/Cantera/src/thermo/IdealGasPhase.h +++ b/Cantera/src/thermo/IdealGasPhase.h @@ -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[ diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index b45d3eb44..73669357e 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -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) diff --git a/Cantera/src/thermo/MetalSHEelectrons.cpp b/Cantera/src/thermo/MetalSHEelectrons.cpp new file mode 100644 index 000000000..221fa5ae9 --- /dev/null +++ b/Cantera/src/thermo/MetalSHEelectrons.cpp @@ -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 +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 + * T and P 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 T and P 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: + * + * + * + * 3.52 + * + * + */ + 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); + } + //==================================================================================================================== + +} diff --git a/Cantera/src/thermo/MetalSHEelectrons.h b/Cantera/src/thermo/MetalSHEelectrons.h new file mode 100644 index 000000000..424e91f2c --- /dev/null +++ b/Cantera/src/thermo/MetalSHEelectrons.h @@ -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 H2 gas at the system pressure + * + * + * Specification of Species Standard %State Properties + * + * 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] + * + * + * Specification of Solution Thermodynamic Properties + * + * All solution properties are obtained from the standard state + * species functions, since there is only one species in the phase. + * + * %Application within %Kinetics Managers + * + * 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. + * + * + * Instanteation of the Class + * + * 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 + * + * + * + * XML Example + * + * 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 + + + + + + + E + + she_electron + + 2.165 + + + + + + + + + E:1 + -1 + + + + 1.172165560E+00, 3.990260375E-03, -9.739075500E-06, 1.007860470E-08, + -3.688058805E-12, -4.589675865E+02, 3.415051190E-01 + + + + + 1.466432895E+00, 4.133039835E-04, -7.320116750E-08, 7.705017950E-12, + -3.444022160E-16, -4.065327985E+02, -5.121644350E-01 + + + + 2.165 + + + +@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 T and P 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 T and P 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 T and P 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 T and P 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 T and P 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 T and P 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 + + + 3.52 + + @endverbatim + * + */ + virtual void setParametersFromXML(const XML_Node& eosdata); + + protected: + + XML_Node *xdef_; + }; + + +} +#endif diff --git a/Cantera/src/thermo/SingleSpeciesTP.cpp b/Cantera/src/thermo/SingleSpeciesTP.cpp index e9980ada2..3bdfe52b3 100644 --- a/Cantera/src/thermo/SingleSpeciesTP.cpp +++ b/Cantera/src/thermo/SingleSpeciesTP.cpp @@ -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; */ diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index c1e33c528..7913d50d4 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -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 diff --git a/Cantera/src/thermo/WaterProps.cpp b/Cantera/src/thermo/WaterProps.cpp index 5501a7129..061499a92 100644 --- a/Cantera/src/thermo/WaterProps.cpp +++ b/Cantera/src/thermo/WaterProps.cpp @@ -495,6 +495,7 @@ namespace Cantera { // Parameters for the viscosityWater() function + //@{ const double H[4] = {1., 0.978197, diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp index d91bce64b..5f8a382de 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.cpp @@ -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; } diff --git a/Cantera/src/thermo/WaterPropsIAPWSphi.h b/Cantera/src/thermo/WaterPropsIAPWSphi.h index 1691214ae..cbd905d21 100644 --- a/Cantera/src/thermo/WaterPropsIAPWSphi.h +++ b/Cantera/src/thermo/WaterPropsIAPWSphi.h @@ -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 diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index adee9aeba..50482d169 100644 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -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;