From ff7e66e5a34e05dff7d678eb255fcf1966090b4e Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 23 Jul 2011 00:30:35 +0000 Subject: [PATCH] This update adds in a Multicomponent Real Gas capability to Cantera. The Object is called RedlichKwongMFTP, and it implements a multicomponent real gas approximation. Most of the functionality has been verified against sample problems involving the CO2-H2O system. --- Cantera/src/base/checkFinite.cpp | 34 +- Cantera/src/numerics/RootFind.h | 1 - Cantera/src/thermo/Makefile.in | 6 +- .../src/thermo/MixedSolventElectrolyte.cpp | 1214 +++++++++++ Cantera/src/thermo/MixedSolventElectrolyte.h | 1025 +++++++++ Cantera/src/thermo/MixtureFugacityTP.cpp | 1367 ++++++++++++ Cantera/src/thermo/MixtureFugacityTP.h | 968 +++++++++ Cantera/src/thermo/NasaThermo.h | 5 +- Cantera/src/thermo/RedlichKwongMFTP.cpp | 1895 +++++++++++++++++ Cantera/src/thermo/RedlichKwongMFTP.h | 843 ++++++++ Cantera/src/thermo/ShomateThermo.h | 5 +- Cantera/src/thermo/SimpleThermo.h | 3 +- Cantera/src/thermo/ThermoFactory.cpp | 11 + Cantera/src/thermo/ThermoPhase.cpp | 4 +- Cantera/src/thermo/ThermoPhase.h | 2 +- Cantera/src/thermo/mix_defs.h | 4 + config.h.in | 3 + configure | 16 +- configure.in | 7 + .../linux.64_sierra_gcc444_python264_numpy | 3 + preconfig | 4 + .../min_python/negATest/negATest_blessed.out | 4 +- 22 files changed, 7390 insertions(+), 34 deletions(-) create mode 100644 Cantera/src/thermo/MixedSolventElectrolyte.cpp create mode 100644 Cantera/src/thermo/MixedSolventElectrolyte.h create mode 100644 Cantera/src/thermo/MixtureFugacityTP.cpp create mode 100644 Cantera/src/thermo/MixtureFugacityTP.h create mode 100644 Cantera/src/thermo/RedlichKwongMFTP.cpp create mode 100644 Cantera/src/thermo/RedlichKwongMFTP.h diff --git a/Cantera/src/base/checkFinite.cpp b/Cantera/src/base/checkFinite.cpp index be4fa224d..443d3d8a7 100644 --- a/Cantera/src/base/checkFinite.cpp +++ b/Cantera/src/base/checkFinite.cpp @@ -50,11 +50,11 @@ namespace mdp { void checkFinite(const double tmp) throw(std::range_error) { if (_finite(tmp)) { if(_isnan(tmp)) { - printf("ERROR: we have encountered a nan!\n"); + printf("checkFinite() ERROR: we have encountered a nan!\n"); } else if (_fpclass(tmp) == _FPCLASS_PINF) { - printf("ERROR: we have encountered a pos inf!\n"); + printf("checkFinite() ERROR: we have encountered a pos inf!\n"); } else { - printf("ERROR: we have encountered a neg inf!\n"); + printf("checkFinite() ERROR: we have encountered a neg inf!\n"); } const std::string s = "checkFinite()"; throw std::range_error(s); @@ -64,11 +64,11 @@ namespace mdp { void checkFinite(const double tmp) throw(std::range_error) { if (! finite(tmp)) { if(isnan(tmp)) { - printf("ERROR: we have encountered a nan!\n"); + printf("checkFinite() ERROR: we have encountered a nan!\n"); } else if (isinf(tmp) == 1) { - printf("ERROR: we have encountered a pos inf!\n"); + printf("checkFinite() ERROR: we have encountered a pos inf!\n"); } else { - printf("ERROR: we have encountered a neg inf!\n"); + printf("checkFinite() ERROR: we have encountered a neg inf!\n"); } const std::string s = "checkFinite()"; throw std::range_error(s); @@ -102,7 +102,7 @@ namespace mdp { checkFinite(tmp); if (fabs(tmp) >= trigger) { char sbuf[64]; - sprintf(sbuf, "checkMagnitude: Trigger %g exceeded: %g\n", trigger, + sprintf(sbuf, "checkMagnitude() ERROR: Trigger %g exceeded: %g\n", trigger, tmp); throw std::range_error(sbuf); } @@ -119,16 +119,16 @@ namespace mdp { void checkZeroFinite(const double tmp) throw(std::range_error) { if ((tmp == 0.0) || (! _finite(tmp))) { if (tmp == 0.0) { - printf("ERROR: we have encountered a zero!\n"); + printf("checkZeroFinite() ERROR: we have encountered a zero!\n"); } else if(_isnan(tmp)) { - printf("ERROR: we have encountered a nan!\n"); + printf("checkZeroFinite() ERROR: we have encountered a nan!\n"); } else if (_fpclass(tmp) == _FPCLASS_PINF) { - printf("ERROR: we have encountered a pos inf!\n"); + printf("checkZeroFinite() ERROR: we have encountered a pos inf!\n"); } else { - printf("ERROR: we have encountered a neg inf!\n"); + printf("checkZeroFinite() ERROR: we have encountered a neg inf!\n"); } char sbuf[64]; - sprintf(sbuf, "checkZeroFinite: zero or indef exceeded: %g\n", + sprintf(sbuf, "checkZeroFinite() ERROR: zero or indef exceeded: %g\n", tmp); throw std::range_error(sbuf); } @@ -137,16 +137,16 @@ void checkZeroFinite(const double tmp) throw(std::range_error) { void checkZeroFinite(const double tmp) throw(std::range_error) { if ((tmp == 0.0) || (! finite(tmp))) { if (tmp == 0.0) { - printf("ERROR: we have encountered a zero!\n"); + printf("checkZeroFinite() ERROR: we have encountered a zero!\n"); } else if(isnan(tmp)) { - printf("ERROR: we have encountered a nan!\n"); + printf("checkZeroFinite() ERROR: we have encountered a nan!\n"); } else if (isinf(tmp) == 1) { - printf("ERROR: we have encountered a pos inf!\n"); + printf("checkZeroFinite() ERROR: we have encountered a pos inf!\n"); } else { - printf("ERROR: we have encountered a neg inf!\n"); + printf("checkZeroFinite() ERROR: we have encountered a neg inf!\n"); } char sbuf[64]; - sprintf(sbuf, "checkZeroFinite: zero or indef exceeded: %g\n", + sprintf(sbuf, "checkZeroFinite() ERROR: zero or indef exceeded: %g\n", tmp); throw std::range_error(sbuf); } diff --git a/Cantera/src/numerics/RootFind.h b/Cantera/src/numerics/RootFind.h index bc5b6e3fc..0392a8799 100644 --- a/Cantera/src/numerics/RootFind.h +++ b/Cantera/src/numerics/RootFind.h @@ -137,7 +137,6 @@ namespace Cantera { */ int solve(doublereal xmin, doublereal xmax, int itmax, doublereal &funcTargetValue, doublereal *xbest); - //! Return the function value /*! * This routine evaluates the following equation. diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index 690877d38..9a933701f 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -40,7 +40,7 @@ THERMO_OBJ = State.o Elements.o Constituents.o Phase.o \ ThermoFactory.o phasereport.o SpeciesThermoInterpType.o \ VPSSMgr.o VPSSMgrFactory.o VPSSMgr_General.o IdealSolnGasVPSS.o \ VPSSMgr_IdealGas.o VPSSMgr_ConstVol.o PDSS_ConstVol.o PDSS_IdealGas.o \ - PDSS_SSVol.o @phase_object_files@ + PDSS_SSVol.o MixtureFugacityTP.o RedlichKwongMFTP.o @phase_object_files@ THERMO_H = State.h Elements.h Constituents.h Phase.h mix_defs.h \ ThermoPhase.h IdealGasPhase.h ConstDensityThermo.h \ @@ -54,7 +54,7 @@ THERMO_H = State.h Elements.h Constituents.h Phase.h mix_defs.h \ EdgePhase.h \ VPSSMgr.h VPSSMgrFactory.h VPSSMgr_General.h IdealSolnGasVPSS.h \ VPSSMgr_IdealGas.h VPSSMgr_ConstVol.h PDSS_ConstVol.h PDSS_IdealGas.h \ - PDSS_SSVol.h @phase_header_files@ + PDSS_SSVol.h MixtureFugacityTP.h RedlichKwongMFTP.h @phase_header_files@ # Extended Cantera Thermodynamics Object Files @@ -91,7 +91,7 @@ CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ) CATHERMO_H = $(THERMO_H) $(ELECTRO_H) $(ISSP_H) -CXX_INCLUDES = -I../base @CXX_INCLUDES@ +CXX_INCLUDES = -I../base -I../numerics @CXX_INCLUDES@ LIB = @buildlib@/libthermo.a DEPENDS = $(CATHERMO_OBJ:.o=.d) diff --git a/Cantera/src/thermo/MixedSolventElectrolyte.cpp b/Cantera/src/thermo/MixedSolventElectrolyte.cpp new file mode 100644 index 000000000..64d4b9bd1 --- /dev/null +++ b/Cantera/src/thermo/MixedSolventElectrolyte.cpp @@ -0,0 +1,1214 @@ +/** + * @file MargulesVPSSTP.cpp + * Definitions for ThermoPhase object for phases which + * employ excess gibbs free energy formulations related to Margules + * expansions (see \ref thermoprops + * and class \link Cantera::MargulesVPSSTP MargulesVPSSTP\endlink). + * + */ +/* + * Copywrite (2009) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Date: 2011-04-14 12:24:13 -0600 (Thu, 14 Apr 2011) $ + * $Revision: 713 $ + */ + + +#include "MixedSolventElectrolyte.h" +#include "ThermoFactory.h" +#include + +using namespace std; + +namespace Cantera { + + static const double xxSmall = 1.0E-150; + /* + * Default constructor. + * + */ + MixedSolventElectrolyte::MixedSolventElectrolyte() : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + } + + /* + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + + */ + MixedSolventElectrolyte::MixedSolventElectrolyte(std::string inputFile, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + constructPhaseFile(inputFile, id); + } + + MixedSolventElectrolyte::MixedSolventElectrolyte(XML_Node& phaseRoot, std::string id) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + constructPhaseXML(phaseRoot, id); + } + + + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor + */ + MixedSolventElectrolyte::MixedSolventElectrolyte(const MixedSolventElectrolyte &b) : + GibbsExcessVPSSTP() + { + MixedSolventElectrolyte::operator=(b); + } + + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + MixedSolventElectrolyte& MixedSolventElectrolyte:: + operator=(const MixedSolventElectrolyte &b) { + if (&b == this) { + return *this; + } + + GibbsExcessVPSSTP::operator=(b); + + numBinaryInteractions_ = b.numBinaryInteractions_ ; + m_HE_b_ij = b.m_HE_b_ij; + m_HE_c_ij = b.m_HE_c_ij; + m_HE_d_ij = b.m_HE_d_ij; + m_SE_b_ij = b.m_SE_b_ij; + m_SE_c_ij = b.m_SE_c_ij; + m_SE_d_ij = b.m_SE_d_ij; + m_VHE_b_ij = b.m_VHE_b_ij; + m_VHE_c_ij = b.m_VHE_c_ij; + m_VHE_d_ij = b.m_VHE_d_ij; + m_VSE_b_ij = b.m_VSE_b_ij; + m_VSE_c_ij = b.m_VSE_c_ij; + m_VSE_d_ij = b.m_VSE_d_ij; + m_pSpecies_A_ij = b.m_pSpecies_A_ij; + m_pSpecies_B_ij = b.m_pSpecies_B_ij; + formMargules_ = b.formMargules_; + formTempModel_ = b.formTempModel_; + + return *this; + } + + /** + * + * ~MixedSolventElectrolyte(): (virtual) + * + * Destructor: does nothing: + * + */ + MixedSolventElectrolyte::~MixedSolventElectrolyte() { + } + + /* + * This routine duplicates the current object and returns + * a pointer to ThermoPhase. + */ + ThermoPhase* + MixedSolventElectrolyte::duplMyselfAsThermoPhase() const { + MixedSolventElectrolyte* mtp = new MixedSolventElectrolyte(*this); + return (ThermoPhase *) mtp; + } + + // Special constructor for a hard-coded problem + /* + * + * LiKCl treating the PseudoBinary layer as passthrough. + * -> test to predict the eutectic and liquidus correctly. + * + */ + MixedSolventElectrolyte::MixedSolventElectrolyte(int testProb) : + GibbsExcessVPSSTP(), + numBinaryInteractions_(0), + formMargules_(0), + formTempModel_(0) + { + + + constructPhaseFile("LiKCl_liquid.xml", ""); + + + numBinaryInteractions_ = 1; + + m_HE_b_ij.resize(1); + m_HE_c_ij.resize(1); + m_HE_d_ij.resize(1); + + m_SE_b_ij.resize(1); + m_SE_c_ij.resize(1); + m_SE_d_ij.resize(1); + + m_VHE_b_ij.resize(1); + m_VHE_c_ij.resize(1); + m_VHE_d_ij.resize(1); + + m_VSE_b_ij.resize(1); + m_VSE_c_ij.resize(1); + m_VSE_d_ij.resize(1); + + m_pSpecies_A_ij.resize(1); + m_pSpecies_B_ij.resize(1); + + + + m_HE_b_ij[0] = -17570E3; + m_HE_c_ij[0] = -377.0E3; + m_HE_d_ij[0] = 0.0; + + m_SE_b_ij[0] = -7.627E3; + m_SE_c_ij[0] = 4.958E3; + m_SE_d_ij[0] = 0.0; + + + int iLiCl = speciesIndex("LiCl(L)"); + if (iLiCl < 0) { + throw CanteraError("MixedSolventElectrolyte test1 constructor", + "Unable to find LiCl(L)"); + } + m_pSpecies_B_ij[0] = iLiCl; + + + int iKCl = speciesIndex("KCl(L)"); + if (iKCl < 0) { + throw CanteraError("MixedSolventElectrolyte test1 constructor", + "Unable to find KCl(L)"); + } + m_pSpecies_A_ij[0] = iKCl; + } + + + /* + * -------------- Utilities ------------------------------- + */ + + + // Equation of state type flag. + /* + * The ThermoPhase base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Known constants defined for this purpose are + * listed in mix_defs.h. The MixedSolventElectrolyte class also returns + * zero, as it is a non-complete class. + */ + int MixedSolventElectrolyte::eosType() const { + return 0; + } + + /* + * Import, construct, and initialize a phase + * specification from an XML tree into the current object. + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void MixedSolventElectrolyte::constructPhaseFile(std::string inputFile, std::string id) { + + if ((int) inputFile.size() == 0) { + throw CanteraError("MixedSolventElectrolyte:constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + std::ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("MixedSolventElectrolyte:constructPhaseFile","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + XML_Node &phaseNode_XML = xml(); + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("MixedSolventElectrolyte:constructPhaseFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + + /* + * Import, construct, and initialize a HMWSoln phase + * specification from an XML tree into the current object. + * + * Most of the work is carried out by the cantera base + * routine, importPhase(). That routine imports all of the + * species and element data, including the standard states + * of the species. + * + * Then, In this routine, we read the information + * particular to the specification of the activity + * coefficient model for the Pitzer parameterization. + * + * We also read information about the molar volumes of the + * standard states if present in the XML file. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void MixedSolventElectrolyte::constructPhaseXML(XML_Node& phaseNode, std::string id) { + string stemp; + if ((int) id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError("MixedSolventElectrolyte::constructPhaseXML", + "phasenode and Id are incompatible"); + } + } + + /* + * Find the Thermo XML node + */ + if (!phaseNode.hasChild("thermo")) { + throw CanteraError("MixedSolventElectrolyte::constructPhaseXML", + "no thermo XML node"); + } + XML_Node& thermoNode = phaseNode.child("thermo"); + + /* + * Make sure that the thermo model is Margules + */ + stemp = thermoNode.attrib("model"); + string formString = lowercase(stemp); + if (formString != "margules") { + throw CanteraError("MixedSolventElectrolyte::constructPhaseXML", + "model name isn't Margules: " + formString); + + } + + /* + * Call the Cantera importPhase() function. This will import + * all of the species into the phase. This will also handle + * all of the solvent and solute standard states + */ + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError("MixedSolventElectrolyte::constructPhaseXML","importPhase failed "); + } + + } + + + /* + * ------------ Molar Thermodynamic Properties ---------------------- + */ + + + /* + * - Activities, Standard States, Activity Concentrations ----------- + */ + + // 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. Note that they may + * or may not have units of concentration --- they might be + * partial pressures, mole fractions, or surface coverages, + * for example. + * + * Here we define the activity concentrations as equal + * to the activities, because the standard concentration is 1. + * + * @param c Output array of generalized concentrations. The + * units depend upon the implementation of the + * reaction rate expressions within the phase. + */ + void MixedSolventElectrolyte::getActivityConcentrations(doublereal* c) const { + getActivities(c); + } + + doublereal MixedSolventElectrolyte::standardConcentration(int k) const { + //err("standardConcentration"); + //return -1.0; + return 1.0; + } + + doublereal MixedSolventElectrolyte::logStandardConc(int k) const { + //err("logStandardConc"); + //return -1.0; + return 0.0; + } + //==================================================================================================================== + // Get the array of non-dimensional molar-based activity coefficients at + // the current solution temperature, pressure, and solution concentration. + /* + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + void MixedSolventElectrolyte::getActivityCoefficients(doublereal* ac) const { + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + + /* + * take the exp of the internally storred coefficients. + */ + for (int k = 0; k < m_kk; k++) { + ac[k] = exp(lnActCoeff_Scaled_[k]); + } + } + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + + + + void MixedSolventElectrolyte::getElectrochemPotentials(doublereal* mu) const { + getChemPotentials(mu); + double ve = Faraday * electricPotential(); + for (int k = 0; k < m_kk; k++) { + mu[k] += ve*charge(k); + } + } + + + void MixedSolventElectrolyte::getChemPotentials(doublereal* mu) const { + doublereal xx; + /* + * First get the standard chemical potentials in + * molar form. + * -> this requires updates of standard state as a function + * of T and P + */ + getStandardChemPotentials(mu); + /* + * Update the activity coefficients + */ + s_update_lnActCoeff(); + /* + * + */ + doublereal RT = GasConstant * temperature(); + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(moleFractions_[k], xxSmall); + mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + } + } + + /// Molar enthalpy. Units: J/kmol. + doublereal MixedSolventElectrolyte::enthalpy_mole() const { + int kk = nSpecies(); + double hbar[kk], h = 0; + getPartialMolarEnthalpies(hbar); + for (int i = 0; i < kk; i++){ + h += moleFractions_[i]*hbar[i]; + } + return h; + } + + /// Molar entropy. Units: J/kmol. + doublereal MixedSolventElectrolyte::entropy_mole() const { + int kk = nSpecies(); + double sbar[kk], s = 0; + getPartialMolarEntropies(sbar); + for (int i = 0; i < kk; i++){ + s += moleFractions_[i]*sbar[i]; + } + return s; + } + + /// Molar heat capacity at constant pressure. Units: J/kmol/K. + doublereal MixedSolventElectrolyte::cp_mole() const { + int kk = nSpecies(); + double cpbar[kk], cp = 0; + getPartialMolarCp(cpbar); + for (int i = 0; i < kk; i++){ + cp += moleFractions_[i]*cpbar[i]; + } + return cp; + } + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + doublereal MixedSolventElectrolyte::cv_mole() const { + return cp_mole() - GasConstant; + } + + // Returns an array of partial molar enthalpies for the species + // in the mixture. + /* + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * molality-based activity coefficent wrt temperature + * + * \f[ + * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + */ + void MixedSolventElectrolyte::getPartialMolarEnthalpies(doublereal* hbar) const { + /* + * Get the nondimensional standard state enthalpies + */ + getEnthalpy_RT(hbar); + /* + * dimensionalize it. + */ + double T = temperature(); + double RT = GasConstant * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] *= RT; + } + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + double RTT = RT * T; + for (int k = 0; k < m_kk; k++) { + hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k]; + } + } + + // Returns an array of partial molar heat capacities for the species + // in the mixture. + /* + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * ??????????? \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + */ + void MixedSolventElectrolyte::getPartialMolarCp(doublereal* cpbar) const { + /* + * Get the nondimensional standard state entropies + */ + getCp_R(cpbar); + double T = temperature(); + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + + for (int k = 0; k < m_kk; k++) { + cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k]; + } + /* + * dimensionalize it. + */ + for (int k = 0; k < m_kk; k++) { + cpbar[k] *= GasConstant; + } + } + + // Returns an array of partial molar entropies for the species + // in the mixture. + /* + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + */ + void MixedSolventElectrolyte::getPartialMolarEntropies(doublereal* sbar) const { + double xx; + /* + * Get the nondimensional standard state entropies + */ + getEntropy_R(sbar); + double T = temperature(); + /* + * Update the activity coefficients, This also update the + * internally storred molalities. + */ + s_update_lnActCoeff(); + s_update_dlnActCoeff_dT(); + + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(moleFractions_[k], xxSmall); + sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k]; + } + /* + * dimensionalize it. + */ + for (int k = 0; k < m_kk; k++) { + sbar[k] *= GasConstant; + } + } + + /* + * ------------ Partial Molar Properties of the Solution ------------ + */ + + // Return an array of partial molar volumes for the + // species in the mixture. Units: m^3/kmol. + /* + * Frequently, for this class of thermodynamics representations, + * the excess Volume due to mixing is zero. Here, we set it as + * a default. It may be overriden in derived classes. + * + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + void MixedSolventElectrolyte::getPartialMolarVolumes(doublereal* vbar) const { + + int iA, iB, iK, delAK, delBK; + double XA, XB, XK, g0 , g1; + double T = temperature(); + + /* + * Get the standard state values in m^3 kmol-1 + */ + getStandardVolumes(vbar); + + + for ( iK = 0; iK < m_kk; iK++ ){ + delAK = 0; + delBK = 0; + XK = moleFractions_[iK]; + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]); + g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]); + + vbar[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1; + } + } + } + + doublereal MixedSolventElectrolyte::err(std::string msg) const { + throw CanteraError("MixedSolventElectrolyte","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 0; + } + + + /* + * @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 MixedSolventElectrolyte::initThermo() { + initLengths(); + GibbsExcessVPSSTP::initThermo(); + } + + + // Initialize lengths of local variables after all species have + // been identified. + void MixedSolventElectrolyte::initLengths() { + m_kk = nSpecies(); + dlnActCoeffdlnN_.resize(m_kk, m_kk); + } + + /* + * initThermoXML() (virtual from ThermoPhase) + * Import and initialize a ThermoPhase object + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, std::string id) { + string subname = "MixedSolventElectrolyte::initThermoXML"; + string stemp; + + /* + * Check on the thermo field. Must have: + * + */ + + XML_Node& thermoNode = phaseNode.child("thermo"); + string mStringa = thermoNode.attrib("model"); + string mString = lowercase(mStringa); + if (mString != "margules") { + throw CanteraError(subname.c_str(), + "Unknown thermo model: " + mStringa); + } + + + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + XML_Node *acNodePtr = 0; + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& acNode = thermoNode.child("activityCoefficients"); + acNodePtr = &acNode; + string mStringa = acNode.attrib("model"); + string mString = lowercase(mStringa); + if (mString != "margules") { + throw CanteraError(subname.c_str(), + "Unknown activity coefficient model: " + mStringa); + } + int n = acNodePtr->nChildren(); + for (int i = 0; i < n; i++) { + XML_Node &xmlACChild = acNodePtr->child(i); + stemp = xmlACChild.name(); + string nodeName = lowercase(stemp); + /* + * Process a binary salt field, or any of the other XML fields + * that make up the Pitzer Database. Entries will be ignored + * if any of the species in the entry isn't in the solution. + */ + if (nodeName == "binaryneutralspeciesparameters") { + readXMLBinarySpecies(xmlACChild); + + } + } + } + + /* + * Go down the chain + */ + GibbsExcessVPSSTP::initThermoXML(phaseNode, id); + + + } + //=================================================================================================================== + + // Update the activity coefficients + /* + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + * + * he = X_A X_B(B + C X_B) + */ + void MixedSolventElectrolyte::s_update_lnActCoeff() const { + int iA, iB, iK, delAK, delBK; + double XA, XB, XK, g0 , g1; + double T = temperature(); + double RT = GasConstant*T; + fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + for (iK = 0; iK < m_kk; iK++) { + XK = moleFractions_[iK]; + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + delAK = 0; + delBK = 0; + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + lnActCoeff_Scaled_[iK] += (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1; + } + } + } + //=================================================================================================================== + // Update the derivative of the log of the activity coefficients wrt T + /* + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + * + * he = X_A X_B(B + C X_B) + */ + void MixedSolventElectrolyte::s_update_dlnActCoeff_dT() const { + int iA, iB, iK, delAK, delBK; + doublereal XA, XB, g0, g1; + doublereal T = temperature(); + doublereal RTT = GasConstant*T*T; + fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); + fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); + for (iK = 0; iK < m_kk; iK++) { + for (int i = 0; i < numBinaryInteractions_; i++) { + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + delAK = 0; + delBK = 0; + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + g0 = -m_HE_b_ij[i] / RTT; + g1 = -m_HE_c_ij[i] / RTT; + double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1; + dlnActCoeffdT_Scaled_[iK] += temp; + d2lnActCoeffdT2_Scaled_[iK] -= 2.0 * temp / T; + } + } + } + //==================================================================================================================== + void MixedSolventElectrolyte::getdlnActCoeffdT(doublereal *dlnActCoeffdT) const { + s_update_dlnActCoeff_dT(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k]; + } + } + //==================================================================================================================== + void MixedSolventElectrolyte::getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const { + s_update_dlnActCoeff_dT(); + for (int k = 0; k < m_kk; k++) { + d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k]; + } + } + //==================================================================================================================== + + // Get the change in activity coefficients w.r.t. change in state (temp, mole fraction, etc.) along + // a line in parameter space or along a line in physical space + /* + * + * @param dTds Input of temperature change along the path + * @param dXds Input vector of changes in mole fraction along the path. length = m_kk + * Along the path length it must be the case that the mole fractions sum to one. + * @param dlnActCoeffds Output vector of the directional derivatives of the + * log Activity Coefficients along the path. length = m_kk + * units are 1/units(s). if s is a physical coordinate then the units are 1/m. + */ + void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const doublereal * const dXds, + doublereal *dlnActCoeffds) const { + + + int iA, iB, iK, delAK, delBK; + double XA, XB, XK, g0 , g1, dXA, dXB; + double T = temperature(); + double RT = GasConstant*T; + + //fvo_zero_dbl_1(dlnActCoeff, m_kk); + s_update_dlnActCoeff_dT(); + + for ( iK = 0; iK < m_kk; iK++ ){ + + XK = moleFractions_[iK]; + dlnActCoeffds[iK] = 0.0; + + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + delAK = 0; + delBK = 0; + + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + dXA = dXds[iA]; + dXB = dXds[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB + + dlnActCoeffdT_Scaled_[iK]*dTds; + } + } + } + //==================================================================================================================== + // Update the derivative of the log of the activity coefficients wrt dlnN + /* + * This function will be called to update the internally stored gradients of the + * logarithm of the activity coefficients. These are used in the determination + * of the diffusion coefficients. + * + * he = X_A X_B(B + C X_B) + */ + void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const { + int iA, iB, iK, delAK, delBK; + double XA, XB, XK, g0 , g1; + double T = temperature(); + double RT = GasConstant*T; + + fvo_zero_dbl_1(dlnActCoeffdlnN_diag_, m_kk); + + for ( iK = 0; iK < m_kk; iK++ ){ + + XK = moleFractions_[iK]; + + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + delAK = 0; + delBK = 0; + + if (iA==iK) delAK = 1; + else if (iB==iK) delBK = 1; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB))); + +// double gfac = g0 + g1 * XB; +// double gggg = (delBK - XB) * g1; + + +// dlnActCoeffdlnN_diag_[iK] += gfac * delAK * ( - XB + delBK); + +// dlnActCoeffdlnN_diag_[iK] += gfac * delBK * ( - XA + delAK); + +// dlnActCoeffdlnN_diag_[iK] += gfac * (2.0 * XA * XB - delAK * XB - XA * delBK); + +// dlnActCoeffdlnN_diag_[iK] += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBK); + +// dlnActCoeffdlnN_diag_[iK] += gggg * ( - 2.0 * XA * XB + delAK * XB + XA * delBK); + +// dlnActCoeffdlnN_diag_[iK] += - g1 * XA * XB * (- XB + delBK); + } + dlnActCoeffdlnN_diag_[iK] = XK*dlnActCoeffdlnN_diag_[iK];//-XK; + } + } + + //==================================================================================================================== + // Update the derivative of the log of the activity coefficients wrt dlnN + /* + * This function will be called to update the internally stored gradients of the + * logarithm of the activity coefficients. These are used in the determination + * of the diffusion coefficients. + * + */ + void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const { + int iA, iB; + doublereal delAK, delBK; + double XA, XB, g0 , g1, XK,XM; + double T = temperature(); + double RT = GasConstant*T; + + doublereal delAM, delBM; + + dlnActCoeffdlnN_.zero(); + + /* + * Loop over the activity coefficient gamma_k + */ + for (int iK = 0; iK < m_kk; iK++) { + XK = moleFractions_[iK]; + for (int iM = 0; iM < m_kk; iM++) { + XM = moleFractions_[iM]; + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + delAK = 0.0; + delBK = 0.0; + delAM = 0.0; + delBM = 0.0; + if (iA==iK) delAK = 1.0; + else if (iB==iK) delBK = 1.0; + if (iA==iM) delAM = 1.0; + else if (iB==iM) delBM = 1.0; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB)); + dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA); + +// double gfac = g0 + g1 * XB; +// double gggg = (delBK - XB) * g1; + + +// dlnActCoeffdlnN_(iK, iM) += gfac * delAK * ( - XB + delBM); + +// dlnActCoeffdlnN_(iK, iM) += gfac * delBK * ( - XA + delAM); + +// dlnActCoeffdlnN_(iK, iM) += gfac * (2.0 * XA * XB - delAM * XB - XA * delBM); + +// dlnActCoeffdlnN_(iK, iM) += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBM); + +// dlnActCoeffdlnN_(iK, iM) += gggg * ( - 2.0 * XA * XB + delAM * XB + XA * delBM); + +// dlnActCoeffdlnN_(iK, iM) += - g1 * XA * XB * (- XB + delBM); + } + dlnActCoeffdlnN_(iK,iM) = XM*dlnActCoeffdlnN_(iK,iM); + } + } + } + //==================================================================================================================== + void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnX_diag() const { + + int iA, iB; + doublereal XA, XB, g0 , g1; + doublereal T = temperature(); + + fvo_zero_dbl_1(dlnActCoeffdlnX_diag_, m_kk); + + doublereal RT = GasConstant * T; + + + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + XA = moleFractions_[iA]; + XB = moleFractions_[iB]; + + g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; + g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + + dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB); + dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB); + } + } + + //==================================================================================================================== + void MixedSolventElectrolyte::getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const { + s_update_dlnActCoeff_dlnN_diag(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k]; + } + } + //==================================================================================================================== + void MixedSolventElectrolyte::getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const { + s_update_dlnActCoeff_dlnX_diag(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k]; + } + } + //==================================================================================================================== + void MixedSolventElectrolyte::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) { + s_update_dlnActCoeff_dlnN(); + double *data = & dlnActCoeffdlnN_(0,0); + for (int k = 0; k < m_kk; k++) { + for (int m = 0; m < m_kk; m++) { + dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m]; + } + } + } + //==================================================================================================================== + void MixedSolventElectrolyte::resizeNumInteractions(const int num) { + numBinaryInteractions_ = num; + m_HE_b_ij.resize(num, 0.0); + m_HE_c_ij.resize(num, 0.0); + m_HE_d_ij.resize(num, 0.0); + m_SE_b_ij.resize(num, 0.0); + m_SE_c_ij.resize(num, 0.0); + m_SE_d_ij.resize(num, 0.0); + m_VHE_b_ij.resize(num, 0.0); + m_VHE_c_ij.resize(num, 0.0); + m_VHE_d_ij.resize(num, 0.0); + m_VSE_b_ij.resize(num, 0.0); + m_VSE_c_ij.resize(num, 0.0); + m_VSE_d_ij.resize(num, 0.0); + + m_pSpecies_A_ij.resize(num, -1); + m_pSpecies_B_ij.resize(num, -1); + + } + //==================================================================================================================== + + /* + * Process an XML node called "binaryNeutralSpeciesParameters" + * This node contains all of the parameters necessary to describe + * the Margules Interaction for a single binary interaction + * This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + */ + void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node &xmLBinarySpecies) { + string xname = xmLBinarySpecies.name(); + if (xname != "binaryNeutralSpeciesParameters") { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", + "Incorrect name for processing this routine: " + xname); + } + double *charge = DATA_PTR(m_speciesCharge); + string stemp; + int nParamsFound; + vector_fp vParams; + string iName = xmLBinarySpecies.attrib("speciesA"); + if (iName == "") { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "no speciesA attrib"); + } + string jName = xmLBinarySpecies.attrib("speciesB"); + if (jName == "") { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "no speciesB attrib"); + } + /* + * Find the index of the species in the current phase. It's not + * an error to not find the species + */ + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + string ispName = speciesName(iSpecies); + if (charge[iSpecies] != 0) { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "speciesA charge problem"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + string jspName = speciesName(jSpecies); + if (charge[jSpecies] != 0) { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "speciesB charge problem"); + } + + resizeNumInteractions(numBinaryInteractions_ + 1); + int iSpot = numBinaryInteractions_ - 1; + m_pSpecies_A_ij[iSpot] = iSpecies; + m_pSpecies_B_ij[iSpot] = jSpecies; + + int num = xmLBinarySpecies.nChildren(); + for (int iChild = 0; iChild < num; iChild++) { + XML_Node &xmlChild = xmLBinarySpecies.child(iChild); + stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + /* + * Process the binary species interaction child elements + */ + if (nodeName == "excessenthalpy") { + /* + * Get the string containing all of the values + */ + ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEnthalpy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_HE_b_ij[iSpot] = vParams[0]; + m_HE_c_ij[iSpot] = vParams[1]; + } + + if (nodeName == "excessentropy") { + /* + * Get the string containing all of the values + */ + ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEntropy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_SE_b_ij[iSpot] = vParams[0]; + m_SE_c_ij[iSpot] = vParams[1]; + } + + if (nodeName == "excessvolume_enthalpy") { + /* + * Get the string containing all of the values + */ + ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_VHE_b_ij[iSpot] = vParams[0]; + m_VHE_c_ij[iSpot] = vParams[1]; + } + + if (nodeName == "excessvolume_entropy") { + /* + * Get the string containing all of the values + */ + ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Entropy for " + ispName + + "::" + jspName, + "wrong number of params found"); + } + m_VSE_b_ij[iSpot] = vParams[0]; + m_VSE_c_ij[iSpot] = vParams[1]; + } + + + } + + } + +} + diff --git a/Cantera/src/thermo/MixedSolventElectrolyte.h b/Cantera/src/thermo/MixedSolventElectrolyte.h new file mode 100644 index 000000000..e91caaca2 --- /dev/null +++ b/Cantera/src/thermo/MixedSolventElectrolyte.h @@ -0,0 +1,1025 @@ +/** + * @file MargulesVPSSTP.h + * Header for intermediate ThermoPhase object for phases which + * employ gibbs excess free energy based formulations + * (see \ref thermoprops + * and class \link Cantera::MargulesVPSSTP MargulesVPSSTP\endlink). + * + * Header file for a derived class of ThermoPhase that handles + * variable pressure standard state methods for calculating + * thermodynamic properties that are further based upon activities + * based on the molality scale. These include most of the methods for + * calculating liquid electrolyte thermodynamics. + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +/* + * $Id: MargulesVPSSTP.h 713 2011-04-14 18:24:13Z hkmoffa $ + */ + +#ifndef CT_MIXEDSOLVENTELECTROLYTEVPSSTP_H +#define CT_MIXEDSOLVENTELECTROLYTEVPSSTP_H + +#include "PseudoBinaryVPSSTP.h" + +namespace Cantera { + + /** + * @ingroup thermoprops + */ + + + //! MixedSolventElectrolyte is a derived class of GibbsExcessVPSSTP that employs + //! the DH and local Marguless approximations for the excess gibbs free energy + /*! + * + * %MargulesVPSSTP derives from class GibbsExcessVPSSTP which is derived + * from VPStandardStateTP, + * and overloads the virtual methods defined there with ones that + * use expressions appropriate for the Margules Excess gibbs free energy + * approximation. + * + * The independent unknowns are pressure, temperature, and mass fraction. + * + * Several concepts are introduced. The first concept is there are temporary + * variables for holding the species standard state values + * of Cp, H, S, G, and V at the + * last temperature and pressure called. These functions are not recalculated + * if a new call is made using the previous temperature and pressure. Currently, + * these variables and the calculation method are handled by the VPSSMgr class, + * for which VPStandardStateTP owns a pointer to. + * + * To support the above functionality, pressure and temperature variables, + * m_plast_ss and m_tlast_ss, are kept which store the last pressure and temperature + * used in the evaluation of standard state properties. + * + * This class is usually used for nearly incompressible phases. For those phases, it + * makes sense to change the equation of state independent variable from + * density to pressure. The variable m_Pcurrent contains the current value of the + * pressure within the phase. + * + * + *
+ *

Specification of Species Standard %State Properties

+ *
+ * + * All species are defined to have standard states that depend upon both + * the temperature and the pressure. The Margules approximation assumes + * symmetric standard states, where all of the standard state assume + * that the species are in pure component states at the temperatue + * and pressure of the solution. I don't think it prevents, however, + * some species from being dilute in the solution. + * + * + *
+ *

Specification of Solution Thermodynamic Properties

+ *
+ * + * The molar excess Gibbs free energy is given by the following formula which is a sum over interactions i. + * Each of the interactions are binary interactions involving two of the species in the phase, denoted, Ai + * and Bi. + * This is the generalization of the Margules formulation for a phase + * that has more than 2 species. + * + * \f[ + * G^E = \sum_i \left( H_{Ei} - T S_{Ei} \right) + * \f] + * \f[ + * H^E_i = n X_{Ai} X_{Bi} \left( h_{o,i} + h_{1,i} X_{Bi} \right) + * \f] + * \f[ + * S^E_i = n X_{Ai} X_{Bi} \left( s_{o,i} + s_{1,i} X_{Bi} \right) + * \f] + * + * where n is the total moles in the solution. + * + * The activity of a species defined in the phase is given by an excess + * Gibbs free energy formulation. + * + * \f[ + * a_k = \gamma_k X_k + * \f] + * + * where + * + * \f[ + * R T \ln( \gamma_k )= \frac{d(n G^E)}{d(n_k)}\Bigg|_{n_i} + * \f] + * + * Taking the derivatives results in the following expression + * + * \f[ + * R T \ln( \gamma_k )= \sum_i \left( \left( \delta_{Ai,k} X_{Bi} + \delta_{Bi,k} X_{Ai} - X_{Ai} X_{Bi} \right) + * \left( g^E_{o,i} + g^E_{1,i} X_{Bi} \right) + + * \left( \delta_{Bi,k} - X_{Bi} \right) X_{Ai} X_{Bi} g^E_{1,i} \right) + * \f] + * where + * \f$ g^E_{o,i} = h_{o,i} - T s_{o,i} \f$ and \f$ g^E_{1,i} = h_{1,i} - T s_{1,i} \f$ + * and where \f$ X_k \f$ is the mole fraction of species k. + * + * This object inherits from the class VPStandardStateTP. Therefore, the specification and + * calculation of all standard state and reference state values are handled at that level. Various functional + * forms for the standard state are permissible. + * The chemical potential for species k is equal to + * + * \f[ + * \mu_k(T,P) = \mu^o_k(T, P) + R T \ln(\gamma_k X_k) + * \f] + * + * The partial molar entropy for species k is given by the following relation, + * + * \f[ + * \tilde{s}_k(T,P) = s^o_k(T,P) - R \ln( \gamma_k X_k ) + * - R T \frac{d \ln(\gamma_k) }{dT} + * \f] + * + * The partial molar enthalpy for species k is given by + * + * \f[ + * \tilde{h}_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + * The partial molar volume for species k is + * + * \f[ + * \tilde V_k(T,P) = V^o_k(T,P) + R T \frac{d \ln(\gamma_k) }{dP} + * \f] + * + * The partial molar Heat Capacity for species k is + * + * \f[ + * \tilde{C}_{p,k}(T,P) = C^o_{p,k}(T,P) - 2 R T \frac{d \ln( \gamma_k )}{dT} + * - R T^2 \frac{d^2 \ln(\gamma_k) }{{dT}^2} + * \f] + * + *
+ *

%Application within %Kinetics Managers

+ *
+ * + * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / + * C^s_k, \f$ where \f$ C^s_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. + * The activity concentration,\f$ C^a_k \f$,is given by the following expression. + * + * \f[ + * C^a_k = C^s_k X_k = \frac{P}{R T} X_k + * \f] + * + * The standard concentration for species k is independent of k and equal to + * + * \f[ + * C^s_k = C^s = \frac{P}{R T} + * \f] + * + * For example, a bulk-phase binary gas reaction between species j and k, producing + * a new gas species l would have the + * following equation for its rate of progress variable, \f$ R^1 \f$, which has + * units of kmol m-3 s-1. + * + * \f[ + * R^1 = k^1 C_j^a C_k^a = k^1 (C^s a_j) (C^s a_k) + * \f] + * where + * \f[ + * C_j^a = C^s a_j \mbox{\quad and \quad} C_k^a = C^s a_k + * \f] + * + * + * \f$ C_j^a \f$ is the activity concentration of species j, and + * \f$ C_k^a \f$ is the activity concentration of species k. \f$ C^s \f$ + * is the standard concentration. \f$ a_j \f$ is + * the activity of species j which is equal to the mole fraction of j. + * + * The reverse rate constant can then be obtained from the law of microscopic reversibility + * and the equilibrium expression for the system. + * + * \f[ + * \frac{a_j a_k}{ a_l} = K_a^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} ) + * \f] + * + * \f$ K_a^{o,1} \f$ is the dimensionless form of the equilibrium constant, associated with + * the pressure dependent standard states \f$ \mu^o_l(T,P) \f$ and their associated activities, + * \f$ a_l \f$, repeated here: + * + * \f[ + * \mu_l(T,P) = \mu^o_l(T, P) + R T \log(a_l) + * \f] + * + * We can switch over to expressing the equilibrium constant in terms of the reference + * state chemical potentials + * + * \f[ + * K_a^{o,1} = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{P} + * \f] + * + * The concentration equilibrium constant, \f$ K_c \f$, may be obtained by changing over + * to activity concentrations. When this is done: + * + * \f[ + * \frac{C^a_j C^a_k}{ C^a_l} = C^o K_a^{o,1} = K_c^1 = + * \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{RT} + * \f] + * + * %Kinetics managers will calculate the concentration equilibrium constant, \f$ K_c \f$, + * using the second and third part of the above expression as a definition for the concentration + * equilibrium constant. + * + * For completeness, the pressure equilibrium constant may be obtained as well + * + * \f[ + * \frac{P_j P_k}{ P_l P_{ref}} = K_p^1 = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) + * \f] + * + * \f$ K_p \f$ is the simplest form of the equilibrium constant for ideal gases. However, it isn't + * necessarily the simplest form of the equilibrium constant for other types of phases; \f$ K_c \f$ is + * used instead because it is completely general. + * + * The reverse rate of progress may be written down as + * \f[ + * R^{-1} = k^{-1} C_l^a = k^{-1} (C^o a_l) + * \f] + * + * 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 + * constant, \f$ K_c \f$. + * + * \f[ + * k^{-1} = k^1 K^1_c + * \f] + * + * \f$k^{-1} \f$ has units of s-1. + * + * + *
+ *

Instantiation of the Class

+ *
+ * + * + * The constructor for this phase is located in the default ThermoFactory + * for %Cantera. A new %IdealGasPhase may be created by the following code + * snippet: + * + * @code + * XML_Node *xc = get_XML_File("silane.xml"); + * XML_Node * const xs = xc->findNameID("phase", "silane"); + * ThermoPhase *silane_tp = newPhase(*xs); + * IdealGasPhase *silaneGas = dynamic_cast (silane_tp); + * @endcode + * + * or by the following constructor: + * + * @code + * XML_Node *xc = get_XML_File("silane.xml"); + * XML_Node * const xs = xc->findNameID("phase", "silane"); + * IdealGasPhase *silaneGas = new IdealGasPhase(*xs); + * @endcode + * + *
+ *

XML Example

+ *
+ * An example of an XML Element named phase setting up a IdealGasPhase + * object named silane is given below. + * + * + * @verbatim + + + Si H He + + H2 H HE SIH4 SI SIH SIH2 SIH3 H3SISIH SI2H6 + H2SISIH2 SI3H8 SI2 SI3 + + + + + + + @endverbatim + * + * The model attribute "IdealGas" of the thermo XML element identifies the phase as + * being of the type handled by the IdealGasPhase object. + * + * @ingroup thermoprops + * + + */ + class MixedSolventElectrolyte : public GibbsExcessVPSSTP { + + public: + + //! Constructor + /*! + * This doesn't do much more than initialize constants with + * default values for water at 25C. Water molecular weight + * comes from the default elements.xml file. It actually + * differs slightly from the IAPWS95 value of 18.015268. However, + * density conservation and therefore element conservation + * is the more important principle to follow. + */ + MixedSolventElectrolyte(); + + //! Construct and initialize a MixedSolventElectrolyte ThermoPhase object + //! directly from an xml input file + /*! + * Working constructors + * + * The two constructors below are the normal way + * the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + * + * @param inputFile Name of the input file containing the phase XML data + * to set up the object + * @param id ID of the phase in the input file. Defaults to the + * empty string. + */ + MixedSolventElectrolyte(std::string inputFile, std::string id = ""); + + //! Construct and initialize a MixedSolventElectrolyte ThermoPhase object + //! directly from an XML database + /*! + * @param phaseRef XML phase node containing the description of the phase + * @param id id attribute containing the name of the phase. + * (default is the empty string) + */ + MixedSolventElectrolyte(XML_Node& phaseRef, std::string id = ""); + + + //! Special constructor for a hard-coded problem + /*! + * + * @param testProb Hard-coded value. Only the value of 1 is + * used. It's for + * a LiKCl system + * -> test to predict the eutectic and liquidus correctly. + */ + MixedSolventElectrolyte(int testProb); + + //! Copy constructor + /*! + * Note this stuff will not work until the underlying phase + * has a working copy constructor + * + * @param b class to be copied + */ + MixedSolventElectrolyte(const MixedSolventElectrolyte& b); + + //! Assignment operator + /*! + * + * @param b class to be copied. + */ + MixedSolventElectrolyte& operator=(const MixedSolventElectrolyte &b); + + //! Destructor + virtual ~MixedSolventElectrolyte(); + + //! Duplication routine for objects which inherit from ThermoPhase. + /*! + * This virtual routine can be used to duplicate thermophase objects + * inherited from ThermoPhase even if the application only has + * a pointer to ThermoPhase to work with. + */ + virtual ThermoPhase *duplMyselfAsThermoPhase() const; + + /** + * + * @name Utilities + * @{ + */ + + + //! Equation of state type flag. + /*! + * The ThermoPhase base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Known constants defined for this purpose are + * listed in mix_defs.h. The MolalityVPSSTP class also returns + * zero, as it is a non-complete class. + */ + virtual int eosType() const; + + //! Initialization of a phase using an xml file + /*! + * This routine is a precursor to + * routine, which does most of the work. + * + * @param inputFile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void constructPhaseFile(std::string inputFile, std::string id); + + //! Import and initialize a phase + //! specification in an XML tree into the current object. + /*! + * Here we read an XML description of the phase. + * We import descriptions of the elements that make up the + * species in a phase. + * We import information about the species, including their + * reference state thermodynamic polynomials. We then freeze + * the state of the species. + * + * Then, we read the species molar volumes from the xml + * tree to finish the initialization. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void constructPhaseXML(XML_Node& phaseNode, std::string id); + + /** + * @} + * @name Molar Thermodynamic Properties + * @{ + */ + + + /** + * @} + * @name Utilities for Solvent ID and Molality + * @{ + */ + + + + + /** + * @} + * @name Mechanical Properties + * @{ + */ + + /** + * @} + * @name Potential Energy + * + * Species may have an additional potential energy due to the + * presence of external gravitation or electric fields. These + * methods allow specifying a potential energy for individual + * species. + * @{ + */ + + /** + * @} + * @name Activities, Standard States, and Activity Concentrations + * + * The activity \f$a_k\f$ of a species in solution is + * related to the chemical potential by \f[ \mu_k = \mu_k^0(T) + * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is + * the chemical potential at unit activity, which depends only + * on temperature and pressure. + * @{ + */ + + //! 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. Note that they may + * or may not have units of concentration --- they might be + * partial pressures, mole fractions, or surface coverages, + * for example. + * + * @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; + + + /** + * The standard concentration \f$ C^0_k \f$ used to normalize + * the generalized concentration. In many cases, this quantity + * will be the same for all species in a phase - for example, + * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this + * reason, this method returns a single value, instead of an + * array. However, for phases in which the standard + * concentration is species-specific (e.g. surface species of + * different sizes), this method may be called with an + * optional parameter indicating the species. + * + * @param k species index. Defaults to zero. + */ + virtual doublereal standardConcentration(int k=0) const; + + /** + * Returns the natural logarithm of the standard + * concentration of the kth species + * + * @param k species index + */ + virtual doublereal logStandardConc(int k=0) const; + + //! Get the array of non-dimensional molar-based activity coefficients at + //! the current solution temperature, pressure, and solution concentration. + /*! + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + virtual void getActivityCoefficients(doublereal* ac) const; + + + + + //@} + /// @name Partial Molar Properties of the Solution + //@{ + + //! Get the species chemical potentials. Units: J/kmol. + /*! + * This function returns a vector of chemical potentials of the + * species in solution at the current temperature, pressure + * and mole fraction of the solution. + * + * @param mu Output vector of species chemical + * potentials. Length: m_kk. Units: J/kmol + */ + virtual void getChemPotentials(doublereal* mu) const; + + /// Molar enthalpy. Units: J/kmol. + virtual doublereal enthalpy_mole() const; + + /// Molar entropy. Units: J/kmol. + virtual doublereal entropy_mole() const; + + /// Molar heat capacity at constant pressure. Units: J/kmol/K. + virtual doublereal cp_mole() const; + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + virtual doublereal cv_mole() const; + + //! Returns an array of partial molar enthalpies for the species + //! in the mixture. + /*! + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * molality-based activity coefficent wrt temperature + * + * \f[ + * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * \f] + * + * @param hbar Vector of returned partial molar enthalpies + * (length m_kk, units = J/kmol) + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + //! Returns an array of partial molar entropies for the species + //! in the mixture. + /*! + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * - R \ln( \gamma_k X_k) + * - R T \frac{d \ln(\gamma_k) }{dT} + * \f] + * + * @param sbar Vector of returned partial molar entropies + * (length m_kk, units = J/kmol/K) + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + //! Returns an array of partial molar entropies for the species + //! in the mixture. + /*! + * Units (J/kmol) + * + * For this phase, the partial molar enthalpies are equal to the + * standard state enthalpies modified by the derivative of the + * activity coefficent wrt temperature + * + * \f[ + * ??????????????? + * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT} + * - R \ln( \gamma_k X_k) + * - R T \frac{d \ln(\gamma_k) }{dT} + * ??????????????? + * \f] + * + * @param cpbar Vector of returned partial molar heat capacities + * (length m_kk, units = J/kmol/K) + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + + //! Return an array of partial molar volumes for the + //! species in the mixture. Units: m^3/kmol. + /*! + * Frequently, for this class of thermodynamics representations, + * the excess Volume due to mixing is zero. Here, we set it as + * a default. It may be overriden in derived classes. + * + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + //! Get the species electrochemical potentials. + /*! + * These are partial molar quantities. + * This method adds a term \f$ Fz_k \phi_k \f$ to the + * to each chemical potential. + * + * Units: J/kmol + * + * @param mu output vector containing the species electrochemical potentials. + * Length: m_kk., units = J/kmol + */ + void getElectrochemPotentials(doublereal* mu) const; + + //! Get the array of temperature second derivatives of the log activity coefficients + /*! + * This function is a virtual class, but it first appears in GibbsExcessVPSSTP + * class and derived classes from GibbsExcessVPSSTP. + * + * units = 1/Kelvin + * + * @param d2lnActCoeffdT2 Output vector of temperature 2nd derivatives of the + * log Activity Coefficients. length = m_kk + * + */ + virtual void getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const; + + //! Get the array of temperature derivatives of the log activity coefficients + /*! + * This function is a virtual class, but it first appears in GibbsExcessVPSSTP + * class and derived classes from GibbsExcessVPSSTP. + * + * units = 1/Kelvin + * + * @param dlnActCoeffdT Output vector of temperature derivatives of the + * log Activity Coefficients. length = m_kk + * + */ + virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) const; + + + + //@} + /// @name Properties of the Standard State of the Species in the Solution + //@{ + + + + //@} + /// @name Thermodynamic Values for the Species Reference States + //@{ + + + /////////////////////////////////////////////////////// + // + // The methods below are not virtual, and should not + // be overloaded. + // + ////////////////////////////////////////////////////// + + /** + * @name Specific Properties + * @{ + */ + + + /** + * @name Setting the State + * + * These methods set all or part of the thermodynamic + * state. + * @{ + */ + + + + //@} + + /** + * @name Chemical Equilibrium + * Routines that implement the Chemical equilibrium capability + * for a single phase, based on the element-potential method. + * @{ + */ + + + + //@} + + + + /// The following methods are used in the process of constructing + /// the phase and setting its parameters from a specification in an + /// input file. They are not normally used in application programs. + /// To see how they are used, see files importCTML.cpp and + /// ThermoFactory.cpp. + + + /*! + * @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(); + + + /** + * Import and initialize a ThermoPhase object + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void initThermoXML(XML_Node& phaseNode, std::string id); + + /** + * @} + * @name Derivatives of Thermodynamic Variables needed for Applications + * @{ + */ + + //! Get the change in activity coefficients w.r.t. change in state (temp, mole fraction, etc.) along + //! a line in parameter space or along a line in physical space + /*! + * + * @param dTds Input of temperature change along the path + * @param dXds Input vector of changes in mole fraction along the path. length = m_kk + * Along the path length it must be the case that the mole fractions sum to one. + * @param dlnActCoeffds Output vector of the directional derivatives of the + * log Activity Coefficients along the path. length = m_kk + * units are 1/units(s). if s is a physical coordinate then the units are 1/m. + */ + virtual void getdlnActCoeffds(const doublereal dTds, const doublereal * const dXds, doublereal *dlnActCoeffds) const; + + //! Get the array of log concentration-like derivatives of the + //! log activity coefficients - diagonal component + /*! + * This function is a virtual method. For ideal mixtures + * (unity activity coefficients), this can return zero. + * Implementations should take the derivative of the + * logarithm of the activity coefficient with respect to the + * logarithm of the mole fraction. + * + * units = dimensionless + * + * @param dlnActCoeffdlnX_diag Output vector of the diagonal component of the log(mole fraction) + * derivatives of the log Activity Coefficients. + * length = m_kk + */ + virtual void getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const; + + //! Get the array of derivatives of the log activity coefficients wrt mole numbers - diagonal only + /*! + * This function is a virtual method. For ideal mixtures + * (unity activity coefficients), this can return zero. + * Implementations should take the derivative of the + * logarithm of the activity coefficient with respect to the + * logarithm of the concentration-like variable (i.e. mole fraction, + * molality, etc.) that represents the standard state. + * + * units = dimensionless + * + * @param dlnActCoeffdlnN_diag Output vector of the diagonal entries for the log(mole fraction) + * derivatives of the log Activity Coefficients. + * length = m_kk + */ + virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const; + + + //! Get the array of derivatives of the log activity coefficients with respect to the ln species mole numbers + /*! + * Implementations should take the derivative of the logarithm of the activity coefficient with respect to a + * log of a species mole number (with all other species mole numbers held constant) + * + * units = 1 / kmol + * + * dlnActCoeffdlnN[ ld * k + m] will contain the derivative of log act_coeff for the mth + * species with respect to the number of moles of the kth species. + * + * \f[ + * \frac{d \ln(\gamma_m) }{d \ln( n_k ) }\Bigg|_{n_i} + * \f] + * + * @param ld Number of rows in the matrix + * @param dlnActCoeffdlnN Output vector of derivatives of the + * log Activity Coefficients. length = m_kk * m_kk + */ + virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) ; + + //@} + + private: + + //! Process an XML node called "binaryNeutralSpeciesParameters" + /*! + * This node contains all of the parameters necessary to describe + * the Margules model for a particular binary interaction. + * This function reads the XML file and writes the coefficients + * it finds to an internal data structures. + * + * @param xmlBinarySpecies Reference to the XML_Node named "binaryNeutralSpeciesParameters" + * containing the binary interaction + */ + void readXMLBinarySpecies(XML_Node &xmlBinarySpecies); + + //! Resize internal arrays within the object that depend upon the number + //! of binary Margules interaction terms + /*! + * @param num Number of binary Margules interaction terms + */ + void resizeNumInteractions(const int num); + + + //! Initialize lengths of local variables after all species have + //! been identified. + void initLengths(); + + //! Update the activity coefficients + /*! + * This function will be called to update the internally storred + * natural logarithm of the activity coefficients + */ + void s_update_lnActCoeff() const; + + //! Update the derivative of the log of the activity coefficients wrt T + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt temperature. + */ + void s_update_dlnActCoeff_dT() const; + + //! Update the derivative of the log of the activity coefficients + //! wrt log(mole fraction) + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the mole fractions. + */ + void s_update_dlnActCoeff_dlnX_diag() const; + + //! Update the derivative of the log of the activity coefficients + //! wrt log(moles) - diagonal only + /*! + * This function will be called to update the internally storred diagonal entries for the + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the moles. + */ + void s_update_dlnActCoeff_dlnN_diag() const; + + //! Update the derivative of the log of the activity coefficients wrt log(moles_m) + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the mole number of species + */ + void s_update_dlnActCoeff_dlnN() const; + + + private: + //! Error function + /*! + * Print an error string and exit + * + * @param msg Message to be printed + */ + doublereal err(std::string msg) const; + + protected: + + + //! number of binary interaction expressions + int numBinaryInteractions_; + + //! Enthalpy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_HE_b_ij; + + //! Enthalpy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_HE_c_ij; + + //! Enthalpy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_HE_d_ij; + + //! Entropy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_SE_b_ij; + + //! Entropy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_SE_c_ij; + + //! Entropy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_SE_d_ij; + + //! Enthalpy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VHE_b_ij; + + //! Enthalpy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VHE_c_ij; + + //! Enthalpy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VHE_d_ij; + + //! Entropy term for the binary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VSE_b_ij; + + //! Entropy term for the ternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VSE_c_ij; + + //! Entropy term for the quaternary mole fraction interaction of the + //! excess gibbs free energy expression + mutable vector_fp m_VSE_d_ij; + + + + //! vector of species indices representing species A in the interaction + /*! + * Each Margules excess Gibbs free energy term involves two species, A and B. + * This vector identifies species A. + */ + vector_int m_pSpecies_A_ij; + + //! vector of species indices representing species B in the interaction + /*! + * Each Margules excess Gibbs free energy term involves two species, A and B. + * This vector identifies species B. + */ + vector_int m_pSpecies_B_ij; + + //! form of the Margules interaction expression + /*! + * Currently there is only one form. + */ + int formMargules_; + + //! form of the temperatuer dependence of the Margules interaction expression + /*! + * Currently there is only one form -> constant wrt temperature. + */ + int formTempModel_; + + + }; + + + +} + +#endif + + + + + diff --git a/Cantera/src/thermo/MixtureFugacityTP.cpp b/Cantera/src/thermo/MixtureFugacityTP.cpp new file mode 100644 index 000000000..21e032fa6 --- /dev/null +++ b/Cantera/src/thermo/MixtureFugacityTP.cpp @@ -0,0 +1,1367 @@ +/** + * @file MixtureFugacityTP.cpp + * Methods file for a derived class of ThermoPhase that handles + * non-ideal mixtures based on the fugacity models (see \ref thermoprops and + * class \link Cantera::MixtureFugacityTP MixtureFugacityTP\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. + */ +/* + * $Author: hkmoffa $ + * $Date: 2010-01-17 12:08:00 -0700 (Sun, 17 Jan 2010) $ + * $Revision: 388 $ + */ + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "MixtureFugacityTP.h" +#include "VPSSMgr.h" +#include "PDSS.h" + + +#ifndef MIN +# define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif +using namespace std; + +namespace Cantera { + + + + //==================================================================================================================== + /* + * Default constructor + */ + MixtureFugacityTP::MixtureFugacityTP() : + ThermoPhase(), + m_Pcurrent(-1.0), + moleFractions_(0), + iState_(FLUID_GAS), + forcedState_(FLUID_UNDEFINED), + m_Tlast_ref(-1.0), + m_logc0(0.0), + m_h0_RT(0), + m_cp0_R(0), + m_g0_RT(0), + m_s0_R(0) + { + } + //==================================================================================================================== + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor. + * + * The copy constructor just calls the assignment operator + * to do the heavy lifting. + */ + MixtureFugacityTP::MixtureFugacityTP(const MixtureFugacityTP &b) : + ThermoPhase(), + m_Pcurrent(-1.0), + moleFractions_(0), + iState_(FLUID_GAS), + forcedState_(FLUID_UNDEFINED), + m_Tlast_ref(-1.0), + m_logc0(0.0), + m_h0_RT(0), + m_cp0_R(0), + m_g0_RT(0), + m_s0_R(0) + { + MixtureFugacityTP::operator=(b); + } + //==================================================================================================================== + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + MixtureFugacityTP& + MixtureFugacityTP::operator=(const MixtureFugacityTP &b) { + if (&b != this) { + /* + * Mostly, this is a passthrough to the underlying + * assignment operator for the ThermoPhase parent object. + */ + ThermoPhase::operator=(b); + /* + * However, we have to handle data that we own. + */ + m_Pcurrent = b.m_Pcurrent; + moleFractions_ = b.moleFractions_; + iState_ = b.iState_; + forcedState_ = b.forcedState_; + m_Tlast_ref = b.m_Tlast_ref; + m_logc0 = b.m_logc0; + m_h0_RT = b.m_h0_RT; + m_cp0_R = b.m_cp0_R; + m_g0_RT = b.m_g0_RT; + m_s0_R = b.m_s0_R; + /* + * The VPSSMgr object contains shallow pointers. Whenever you have shallow + * pointers, they have to be fixed up to point to the correct objects refering + * back to this ThermoPhase's properties. + */ + //m_VPSS_ptr->initAllPtrs(this, m_spthermo); + /* + * The PDSS objects contains shallow pointers. Whenever you have shallow + * pointers, they have to be fixed up to point to the correct objects refering + * back to this ThermoPhase's properties. This function also sets m_VPSS_ptr + * so it occurs after m_VPSS_ptr is set. + */ + + /* + * Ok, the VPSSMgr object is ready for business. + * We need to resync the temperature and the pressure of the new standard states + * with what is storred in this object. + */ + // m_VPSS_ptr->setState_TP(m_Tlast_ss, m_Plast_ss); + } + return *this; + } + //==================================================================================================================== + /* + * ~MixtureFugacityTP(): (virtual) + * + */ + MixtureFugacityTP::~MixtureFugacityTP() { + + } + + /* + * Duplication function. + * This calls the copy constructor for this object. + */ + ThermoPhase* MixtureFugacityTP::duplMyselfAsThermoPhase() const { + MixtureFugacityTP* vptp = new MixtureFugacityTP(*this); + return (ThermoPhase *) vptp; + } + //==================================================================================================================== + // This method returns the convention used in specification + // of the standard state, of which there are currently two, + // temperature based, and variable pressure based. + /* + * Currently, there are two standard state conventions: + * - Temperature-based activities + * cSS_CONVENTION_TEMPERATURE 0 + * - default + * + * - Variable Pressure and Temperature -based activities + * cSS_CONVENTION_VPSS 1 + */ + int MixtureFugacityTP::standardStateConvention() const { + return cSS_CONVENTION_TEMPERATURE; + } + //==================================================================================================================== + // Set the solution branch to force the ThermoPhase to exist on one branch or another + /* + * @param solnBranch Branch that the solution is restricted to. + * the value -1 means gas. The value -2 means unrestricted. + * Values of zero or greater refer to species dominated condensed phases. + */ + void MixtureFugacityTP::setForcedSolutionBranch(int solnBranch) { + forcedState_ = solnBranch; + } + //==================================================================================================================== + // Report the solution branch which the solution is restricted to + /* + * @return Branch that the solution is restricted to. + * the value -1 means gas. The value -2 means unrestricted. + * Values of zero or greater refer to species dominated condensed phases. + */ + int MixtureFugacityTP::forcedSolutionBranch() const { + return forcedState_; + } + //==================================================================================================================== + // Report the solution branch which the solution is actually on + /* + * @return Branch that the solution is restricted to. + * the value -1 means gas. The value -2 means superfluid.. + * Values of zero or greater refer to species dominated condensed phases. + */ + int MixtureFugacityTP::reportSolnBranchActual() const { + return iState_; + } + //==================================================================================================================== + + /* + * ------------Molar Thermodynamic Properties ------------------------- + */ + //==================================================================================================================== + + doublereal MixtureFugacityTP::err(std::string msg) const { + throw CanteraError("MixtureFugacityTP","Base class method " + +msg+" called. Equation of state type: "+int2str(eosType())); + return 0; + } + //==================================================================================================================== + /* + * ---- Partial Molar Properties of the Solution ----------------- + */ + //==================================================================================================================== + /* + * Get the array of non-dimensional species chemical potentials + * These are partial molar Gibbs free energies. + * \f$ \mu_k / \hat R T \f$. + * Units: unitless + * + * We close the loop on this function, here, calling + * getChemPotentials() and then dividing by RT. + */ + void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const{ + getChemPotentials(muRT); + doublereal invRT = 1.0 / _RT(); + for (int k = 0; k < m_kk; k++) { + muRT[k] *= invRT; + } + } + //==================================================================================================================== + /* + * ----- Thermodynamic Values for the Species Standard States States ---- + */ + void MixtureFugacityTP::getStandardChemPotentials(doublereal* g) const { + _updateReferenceStateThermo(); + copy(m_g0_RT.begin(), m_g0_RT.end(), g); + doublereal RT = _RT(); + double tmp = log (pressure() /m_spthermo->refPressure()); + for (int k = 0; k < m_kk; k++) { + g[k] = RT * (g[k] + tmp); + } + } + //==================================================================================================================== + void MixtureFugacityTP::getEnthalpy_RT(doublereal* hrt) const { + getEnthalpy_RT_ref(hrt); + } + //================================================================================================ +#ifdef H298MODIFY_CAPABILITY + // Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1) + /* + * The 298K heat of formation is defined as the enthalpy change to create the standard state + * of the species from its constituent elements in their standard states at 298 K and 1 bar. + * + * @param k Species k + * @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar + */ + void MixtureFugacityTP::modifyOneHf298SS(const int k, const doublereal Hf298New) { + m_spthermo->modifyOneHf298(k, Hf298New); + m_Tlast_ref += 0.0001234; + } +#endif + //==================================================================================================================== + /* + * Get the array of nondimensional entropy functions for the + * standard state species + * at the current T and P of the solution. + */ + void MixtureFugacityTP::getEntropy_R(doublereal* sr) const { + _updateReferenceStateThermo(); + copy(m_s0_R.begin(), m_s0_R.end(), sr); + double tmp = log (pressure() /m_spthermo->refPressure()); + for (int k = 0; k < m_kk; k++) { + sr[k] -= tmp; + } + } + //==================================================================================================================== + /* + * Get the nondimensional gibbs function for the species + * standard states at the current T and P of the solution. + */ + void MixtureFugacityTP::getGibbs_RT(doublereal* grt) const { + _updateReferenceStateThermo(); + copy(m_g0_RT.begin(), m_g0_RT.end(), grt); + double tmp = log (pressure() /m_spthermo->refPressure()); + for (int k = 0; k < m_kk; k++) { + grt[k] += tmp; + } + } + //==================================================================================================================== + /* + * get the pure Gibbs free energies of each species assuming + * it is in its standard state. This is the same as + * getStandardChemPotentials(). + */ + void MixtureFugacityTP::getPureGibbs(doublereal* g) const { + _updateReferenceStateThermo(); + scale(m_g0_RT.begin(), m_g0_RT.end(), g, _RT()); + double tmp = log (pressure() /m_spthermo->refPressure()); + tmp *= _RT(); + for (int k = 0; k < m_kk; k++) { + g[k] += tmp; + } + } + //==================================================================================================================== + /* + * Returns the vector of nondimensional + * internal Energies of the standard state at the current temperature + * and pressure of the solution for each species. + */ + void MixtureFugacityTP::getIntEnergy_RT(doublereal* urt) const { + _updateReferenceStateThermo(); + copy(m_h0_RT.begin(), m_h0_RT.end(), urt); + doublereal p = pressure(); + doublereal tmp = p / _RT(); + doublereal v0 = _RT() / p; + for (int i = 0; i < m_kk; i++) { + urt[i] -= tmp * v0; + } + } + //==================================================================================================================== + /* + * Get the nondimensional heat capacity at constant pressure + * function for the species + * standard states at the current T and P of the solution. + */ + void MixtureFugacityTP::getCp_R(doublereal* cpr) const { + _updateReferenceStateThermo(); + copy(m_cp0_R.begin(), m_cp0_R.end(), cpr); + } + //==================================================================================================================== + /* + * Get the molar volumes of the species standard states at the current + * T and P of the solution. + * units = m^3 / kmol + * + * @param vol Output vector containing the standard state volumes. + * Length: m_kk. + */ + void MixtureFugacityTP::getStandardVolumes(doublereal *vol) const { + _updateReferenceStateThermo(); + doublereal v0 = _RT() / pressure(); + for (int i = 0; i < m_kk; i++) { + vol[i]= v0; + } + } + //==================================================================================================================== + /* + * ----- Thermodynamic Values for the Species Reference States ---- + */ + + /* + * Returns the vector of nondimensional enthalpies of the + * reference state at the current temperature of the solution and + * the reference pressure for the species. + */ + void MixtureFugacityTP::getEnthalpy_RT_ref(doublereal *hrt) const { + _updateReferenceStateThermo(); + copy(m_h0_RT.begin(), m_h0_RT.end(), hrt); + } + //==================================================================================================================== + /* + * Returns the vector of nondimensional + * enthalpies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + */ + void MixtureFugacityTP::getGibbs_RT_ref(doublereal *grt) const { + _updateReferenceStateThermo(); + copy(m_g0_RT.begin(), m_g0_RT.end(), grt); + } + //==================================================================================================================== + /* + * Returns the vector of the + * gibbs function of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * units = J/kmol + * + * This is filled in here so that derived classes don't have to + * take care of it. + */ + void MixtureFugacityTP::getGibbs_ref(doublereal *g) const { + const array_fp& gibbsrt = gibbs_RT_ref(); + scale(gibbsrt.begin(), gibbsrt.end(), g, _RT()); + } + //==================================================================================================================== + const vector_fp & MixtureFugacityTP::gibbs_RT_ref() const { + _updateReferenceStateThermo(); + return m_g0_RT; + } + //==================================================================================================================== + /* + * Returns the vector of nondimensional + * entropies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + */ + void MixtureFugacityTP::getEntropy_R_ref(doublereal *er) const { + _updateReferenceStateThermo(); + copy(m_s0_R.begin(), m_s0_R.end(), er); + return; + } + //==================================================================================================================== + /* + * Returns the vector of nondimensional + * constant pressure heat capacities of the reference state + * at the current temperature of the solution + * and reference pressure for the species. + */ + void MixtureFugacityTP::getCp_R_ref(doublereal *cpr) const { + _updateReferenceStateThermo(); + copy(m_cp0_R.begin(), m_cp0_R.end(), cpr); + } + //==================================================================================================================== + /* + * Get the molar volumes of the species reference states at the current + * T and reference pressure of the solution. + * + * units = m^3 / kmol + */ + void MixtureFugacityTP::getStandardVolumes_ref(doublereal *vol) const { + _updateReferenceStateThermo(); + double pp = refPressure(); + doublereal v0 = _RT() / pp; + for (int i = 0; i < m_kk; i++) { + vol[i]= v0; + } + } + //==================================================================================================================== + // Set the initial state of the phase to the conditions specified in the state XML element. + /* + * + * This method sets the temperature, pressure, and mole fraction vector to a set default value. + * We modify the default behavior here so that TP is evaluated at the same time. + * + * @param state AN XML_Node object corresponding to + * the "state" entry for this phase in the + * input file. + */ + void MixtureFugacityTP::setStateFromXML(const XML_Node& state) { + int doTP = 0; + string comp = ctml::getChildValue(state,"moleFractions"); + if (comp != "") { + // not overloaded in current object -> phase state is not calculated. + setMoleFractionsByName(comp); + doTP = 1; + } else { + comp = ctml::getChildValue(state,"massFractions"); + if (comp != "") { + // not overloaded in current object -> phase state is not calculated. + setMassFractionsByName(comp); + doTP = 1; + } + } + double t = temperature(); + if (state.hasChild("temperature")) { + t = ctml::getFloat(state, "temperature", "temperature"); + doTP = 1; + } + if (state.hasChild("pressure")) { + double p = ctml::getFloat(state, "pressure", "pressure"); + setState_TP(t, p); + } else if (state.hasChild("density")) { + double rho = ctml::getFloat(state, "density", "density"); + setState_TR(t, rho); + } else if (doTP) { + double rho = State::density(); + setState_TR(t, rho); + } + } + //==================================================================================================================== + /* + * Perform initializations after all species have been + * added. + */ + void MixtureFugacityTP::initThermo() { + initLengths(); + ThermoPhase::initThermo(); + } + //==================================================================================================================== + /* + * Initialize the internal lengths. + * (this is not a virtual function) + */ + void MixtureFugacityTP::initLengths() { + m_kk = nSpecies(); + moleFractions_.resize(m_kk, 0.0); + moleFractions_[0] = 1.0; + m_h0_RT.resize(m_kk, 0.0); + m_cp0_R.resize(m_kk, 0.0); + m_g0_RT.resize(m_kk, 0.0); + m_s0_R.resize(m_kk, 0.0); + } + //==================================================================================================================== + void MixtureFugacityTP::setTemperature(const doublereal temp) { + _updateReferenceStateThermo(); + setState_TR(temperature(), density()); + } + //==================================================================================================================== + void MixtureFugacityTP::setPressure(doublereal p) { + setState_TP(temperature(), p); + // double chemPot[5], mf[5]; + // getMoleFractions(mf); + // getChemPotentials(chemPot); + // for (int i = 0; i < m_kk; i++) { + // printf(" MixFug:setPres: mu(%d = %g) = %18.8g\n", i, mf[i], chemPot[i]); + // } + } + //==================================================================================================================== + void MixtureFugacityTP::setMassFractions(const doublereal* const y) { + State::setMassFractions(y); + getMoleFractions(DATA_PTR(moleFractions_)); + } + //==================================================================================================================== + void MixtureFugacityTP::setMassFractions_NoNorm(const doublereal* const y) { + State::setMassFractions_NoNorm(y); + getMoleFractions(DATA_PTR(moleFractions_)); + } + //==================================================================================================================== + void MixtureFugacityTP::setMoleFractions(const doublereal* const x) { + State::setMoleFractions(x); + getMoleFractions(DATA_PTR(moleFractions_)); + } + //==================================================================================================================== + void MixtureFugacityTP::setMoleFractions_NoNorm(const doublereal* const x) { + State::setMoleFractions_NoNorm(x); + getMoleFractions(DATA_PTR(moleFractions_)); + } + //==================================================================================================================== + void MixtureFugacityTP::setConcentrations(const doublereal* const c) { + State::setConcentrations(c); + getMoleFractions(DATA_PTR(moleFractions_)); + } + //==================================================================================================================== + void MixtureFugacityTP::setMoleFractions_NoState(const doublereal* const x) { + State::setMoleFractions(x); + getMoleFractions(DATA_PTR(moleFractions_)); + updateMixingExpressions(); + } + //==================================================================================================================== + void MixtureFugacityTP::calcDensity() { + err("MixtureFugacityTP::calcDensity() called, but EOS for phase is not known"); + } + //==================================================================================================================== + + void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) { + /* + * A pretty tricky algorithm is needed here, due to problems involving + * standard states of real fluids. For those cases you need + * to combine the T and P specification for the standard state, or else + * you may venture into the forbidden zone, especially when nearing the + * triple point. + * Therefore, we need to do the standard state thermo calc with the + * (t, pres) combo. + */ + getMoleFractions(DATA_PTR(moleFractions_)); + + + State::setTemperature(t); + _updateReferenceStateThermo(); + // Depends on the mole fractions and the temperature + updateMixingExpressions(); + // setPressure(pres); + m_Pcurrent = pres; + // double mmw = meanMolecularWeight(); + + if (forcedState_ == FLUID_UNDEFINED) { + double rhoNow = State::density(); + double rho = densityCalc(t, pres, iState_, rhoNow); + if (rho > 0.0) { + State::setDensity(rho); + m_Pcurrent = pres; + iState_ = phaseState(true); + } else { + if (rho < -1.5) { + rho = densityCalc(t, pres, FLUID_UNDEFINED , rhoNow); + if (rho > 0.0) { + State::setDensity(rho); + m_Pcurrent = pres; + iState_ = phaseState(true); + } else { + throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho"); + } + } else { + throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho"); + } + } + + + + } else if (forcedState_ == FLUID_GAS) { + // Normal density calculation + if (iState_ < FLUID_LIQUID_0) { + double rhoNow = State::density(); + double rho = densityCalc(t, pres, iState_, rhoNow); + if (rho > 0.0) { + State::setDensity(rho); + m_Pcurrent = pres; + iState_ = phaseState(true); + if (iState_ >= FLUID_LIQUID_0) { + throw CanteraError("MixtureFugacityTP::setState_TP()", "wrong state"); + } + } else { + throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho"); + } + + } + + + } else if (forcedState_ > FLUID_LIQUID_0) { + if (iState_ >= FLUID_LIQUID_0) { + double rhoNow = State::density(); + double rho = densityCalc(t, pres, iState_, rhoNow); + if (rho > 0.0) { + State::setDensity(rho); + m_Pcurrent = pres; + iState_ = phaseState(true); + if (iState_ == FLUID_GAS) { + throw CanteraError("MixtureFugacityTP::setState_TP()", "wrong state"); + } + } else { + throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho"); + } + + } + } + + + + //setTemperature(t); + //setPressure(pres); + //calcDensity(); + } + //==================================================================================================================== + // Set the internally storred temperature (K) and density (kg/m^3) + /* + * This overrides the default behavior. In addition to just storring the state in the object, we need to do + * an equation of state calculation and figure out what phase state we are in. + * + * @param t Temperature in kelvin + * @param rho Density (kg/m^3) + */ + void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho) { + getMoleFractions(DATA_PTR(moleFractions_)); + State::setTemperature(T); + _updateReferenceStateThermo(); + State::setDensity(rho); + doublereal mv = molarVolume(); + // depends on mole fraction and temperature + updateMixingExpressions(); + + m_Pcurrent = pressureCalc(T, mv); + iState_ = phaseState(true); + + // printf("setState_TR: state at T = %g, rho = %g, mv = %g, P = %20.13g, iState = %d\n", T, rho, mv, m_Pcurrent, iState_); + } + + //==================================================================================================================== + // Set the temperature (K), pressure (Pa), and mole fractions. + /* + * Note, the mole fractions are set first before the pressure is set. + * Setting the pressure may involve the solution of a nonlinear equation. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + * @param x Vector of mole fractions. + * Length is equal to m_kk. + */ + void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublereal* x) { + setMoleFractions_NoState(x); + setState_TP(t,p); + } + //==================================================================================================================== + /* + * Import and initialize a ThermoPhase object + * + * param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + * + * This routine initializes the lengths in the current object and + * then calls the parent routine. + */ + void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, std::string id) { + MixtureFugacityTP::initLengths(); + + //m_VPSS_ptr->initThermo(); + + // m_VPSS_ptr->initThermoXML(phaseNode, id); + ThermoPhase::initThermoXML(phaseNode, id); + } + //==================================================================================================================== + doublereal MixtureFugacityTP::z() const { + doublereal p = pressure(); + doublereal rho = density(); + doublereal mmw = meanMolecularWeight(); + doublereal molarV = mmw / rho; + doublereal rt = _RT(); + doublereal zz = p * molarV / rt; + return zz; + } + //==================================================================================================================== + doublereal MixtureFugacityTP::sresid() const { + throw CanteraError("MixtureFugacityTP::sresid()", "Base Class: not implemented"); + return 0.0; + } + //==================================================================================================================== + doublereal MixtureFugacityTP::hresid() const { + throw CanteraError("MixtureFugacityTP::hresid()", "Base Class: not implemented"); + return 0.0; + } + //==================================================================================================================== + doublereal MixtureFugacityTP::psatEst(doublereal TKelvin) const { + doublereal tcrit = critTemperature(); + doublereal pcrit = critPressure(); + doublereal tt = tcrit/TKelvin; + if (tt < 1.0) { + return pcrit; + } + doublereal lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918; + return pcrit*exp(lpr); + } + //==================================================================================================================== + doublereal MixtureFugacityTP::liquidVolEst(doublereal TKelvin, doublereal &pres) const { + throw CanteraError("MixtureFugacityTP::liquidVolEst()", "unimplemented"); + return 0.0; + } + //==================================================================================================================== + /* + * Calculates the density given the temperature and the pressure, + * and a guess at the density. Note, below T_c, this is a + * multivalued function. This function assumes that the phase is on one side of the vapor dome + * or the other. It does not allow for crosses of the vapor dome. + * + * parameters: + * temperature: Kelvin + * pressure : Pressure in Pascals (Newton/m**2) + * phase : guessed phase of water + * : -1: no guessed phase + * rhoguess : guessed density of the water + * + * -1.0 no guessed density + * + * If a problem is encountered, a negative 1 is returned. + * + * @TODO make this a const function + */ + doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa, + int phase, doublereal rhoguess) { + double tcrit = critTemperature(); + doublereal mmw = meanMolecularWeight(); + // double pcrit = critPressure(); + // doublereal deltaGuess = 0.0; + if (rhoguess == -1.0) { + if (phase != -1) { + if (TKelvin > tcrit) { + rhoguess = presPa * mmw / (GasConstant * TKelvin); + } else { + if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) { + rhoguess = presPa * mmw / (GasConstant * TKelvin); + } else if (phase >= FLUID_LIQUID_0) { + double lqvol = liquidVolEst(TKelvin, presPa); + rhoguess = mmw / lqvol; + } + } + } else { + /* + * Assume the Gas phase initial guess, if nothing is + * specified to the routine + */ + rhoguess = presPa * mmw / (GasConstant * TKelvin); + } + + } + + double molarVolBase = mmw / rhoguess; + double molarVolLast = molarVolBase; + double vc = mmw / critDensity(); + /* + * molar volume of the spinodal at the current temperature and mole fractions. this will + * be updated as we go. + */ + double molarVolSpinodal = vc; + doublereal pcheck = 1.0E-30 + 1.0E-8 * presPa; + doublereal presBase, dpdVBase, delMV; + bool conv = false; + /* + * We start on one side of the vc and stick with that side + */ + bool gasSide = molarVolBase > vc; + if (gasSide) { + molarVolLast = (GasConstant * TKelvin)/presPa; + } else { + molarVolLast = liquidVolEst(TKelvin, presPa); + } + + /* + * OK, now we do a small solve to calculate the molar volume given the T,P value. + * The algorithm is taken from dfind() + */ + for (int n = 0; n < 200; n++) { + + /* + * Calculate the predicted reduced pressure, pred0, based on the + * current tau and dd. + */ + + /* + * Calculate the derivative of the predicted pressure + * wrt the molar volume. + * This routine also returns the pressure, presBase + */ + dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase); + + /* + * If dpdV is positve, then we are in the middle of the + * 2 phase region and beyond the spinodal stability curve. We need to adjust + * the initial guess outwards and start a new iteration. + */ + if (dpdVBase >= 0.0) { + if (TKelvin > tcrit) { + throw CanteraError("", "confused"); + } + /* + * TODO Spawn a calculation for the value of the spinodal point that is + * very accurate. Answer the question as to wethera solution is + * possible on the current side of the vapor dome. + */ + if (gasSide) { + if (molarVolBase >= vc) { + molarVolSpinodal = molarVolBase; + molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal); + } else { + molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal); + } + } else { + if (molarVolBase <= vc) { + molarVolSpinodal = molarVolBase; + molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal); + } else { + molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal); + } + } + continue; + } + + /* + * Check for convergence + */ + if (fabs(presBase-presPa) < pcheck) { + conv = true; + break; + } + + /* + * Dampen and crop the update + */ + doublereal dpdV = dpdVBase; + if (n < 10) { + dpdV = dpdVBase * 1.5; + } + // if (dpdV > -0.001) dpdV = -0.001; + + /* + * Formulate the update to the molar volume by + * Newton's method. Then, crop it to a max value + * of 0.1 times the current volume + */ + delMV = - (presBase - presPa) / dpdV; + if (!gasSide || delMV < 0.0) { + if (fabs(delMV) > 0.2 * molarVolBase) { + delMV = delMV / fabs(delMV) * 0.2 * molarVolBase; + } + } + /* + * Only go 1/10 the way towards the spinodal at any one time. + */ + if (TKelvin < tcrit) { + if (gasSide) { + if (delMV < 0.0) { + if (-delMV > 0.5 * (molarVolBase - molarVolSpinodal)) { + delMV = - 0.5 * (molarVolBase - molarVolSpinodal); + } + } + } else { + if (delMV > 0.0) { + if (delMV > 0.5 * (molarVolSpinodal - molarVolBase)) { + delMV = 0.5 * (molarVolSpinodal - molarVolBase); + } + } + } + } + /* + * updated the molar volume value + */ + molarVolLast = molarVolBase; + molarVolBase += delMV; + + + if (fabs(delMV/molarVolBase) < 1.0E-14) { + conv = true; + break; + } + + /* + * Check for negative molar volumes + */ + if (molarVolBase <= 0.0) { + molarVolBase = MIN(1.0E-30, fabs(delMV*1.0E-4)); + } + + } + + + /* + * Check for convergence, and return 0.0 if it wasn't achieved. + */ + double densBase = 0.0; + if (! conv) { + molarVolBase = 0.0; + throw CanteraError("MixtureFugacityTP::densityCalc()", "Process didnot converge"); + } else { + densBase = mmw / molarVolBase; + } + return densBase; + } + //==================================================================================================================== + void MixtureFugacityTP::updateMixingExpressions() { + + } + //==================================================================================================================== + MixtureFugacityTP::spinodalFunc::spinodalFunc(MixtureFugacityTP *tp) : + ResidEval(), + m_tp(tp) + { + } + //==================================================================================================================== + int MixtureFugacityTP::spinodalFunc::evalSS(const doublereal t, const doublereal * const y, + doublereal * const r) { + int status = 0; + doublereal molarVol = y[0]; + doublereal tt = m_tp->temperature(); + doublereal pp; + doublereal val = m_tp->dpdVCalc(tt, molarVol, pp); + r[0] = val; + return status; + } + //==================================================================================================================== + // Utility routine in the calculation of the saturation pressure + /* + * Private routine + * + * @param TKelvin temperature (kelvin) + * @param pres pressure (Pascal) + * @param densLiq Output density of liquid + * @param densGas output density of gas + * @param delGRT output delGRT + * + * @return Returns zero if both the gas and the liquid states are found for a given pressure. + + */ + int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal &densLiqGuess, + doublereal &densGasGuess, doublereal &liqGRT, doublereal &gasGRT) { + + int retn = 0; + doublereal densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess); + if (densLiq <= 0.0) { + // throw Cantera::CanteraError("MixtureFugacityTP::corr0", + // "Error occurred trying to find liquid density at (T,P) = " + // + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres)); + retn = -1; + } else { + densLiqGuess = densLiq; + setState_TR(TKelvin, densLiq); + liqGRT = gibbs_mole() / _RT(); + } + + doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess); + if (densGas <= 0.0) { + //throw Cantera::CanteraError("MixtureFugacityTP::corr0", + // "Error occurred trying to find gas density at (T,P) = " + // + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres)); + if (retn == -1) { + throw Cantera::CanteraError("MixtureFugacityTP::corr0", + "Error occurred trying to find gas density at (T,P) = " + + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres)); + } + retn = -2; + } else { + densGasGuess = densGas; + setState_TR(TKelvin, densGas); + gasGRT = gibbs_mole() / _RT(); + } + // delGRT = gibbsLiqRT - gibbsGasRT; + return retn; + } + //==================================================================================================================== + // Returns the Phase State flag for the current state of the object + /* + * @param checkState If true, this function does a complete check to see where + * in paramters space we are + * + * There are three values: + * WATER_GAS below the critical temperature but below the critical density + * WATER_LIQUID below the critical temperature but above the critical density + * WATER_SUPERCRIT above the critical temperature + */ + int MixtureFugacityTP::phaseState(bool checkState) const { + int state; + if (checkState) { + double t = temperature(); + double tcrit = critTemperature(); + double rhocrit = critDensity(); + if (t >= tcrit) { + state = FLUID_SUPERCRIT; + return state; + } + double tmid = tcrit - 100.; + if (tmid < 0.0) { + tmid = tcrit / 2.0; + } + double pp = psatEst(tmid); + double mmw = meanMolecularWeight(); + double molVolLiqTmid = liquidVolEst(tmid, pp); + double molVolGasTmid = GasConstant * tmid / (pp); + double densLiqTmid = mmw / molVolLiqTmid; + double densGasTmid = mmw / molVolGasTmid; + double densMidTmid = 0.5 * (densLiqTmid + densGasTmid); + doublereal rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid); + + double rho = density(); + int iStateGuess = FLUID_LIQUID_0; + if (rho < rhoMid) { + iStateGuess = FLUID_GAS; + } + double molarVol = mmw / rho; + double presCalc; + + double dpdv = dpdVCalc(t, molarVol, presCalc); + if (dpdv < 0.0) { + state = iStateGuess; + } else { + state = FLUID_UNSTABLE; + } + + } + return state; + } + //==================================================================================================================== + // Return the value of the density at the liquid spinodal point (on the liquid side) + // for the current temperature. + /* + * @return returns the density with units of kg m-3 + */ + doublereal MixtureFugacityTP::densSpinodalLiquid() const { + throw CanteraError("", "unimplmented"); + return 0.0; + } + //==================================================================================================================== + // Return the value of the density at the gas spinodal point (on the gas side) + // for the current temperature. + /* + * @return returns the density with units of kg m-3 + */ + doublereal MixtureFugacityTP::densSpinodalGas() const { + throw CanteraError("", "unimplmented"); + return 0.0; + } + //==================================================================================================================== + // Calculate the saturation pressure at the current mixture content for the given temperature + /* + * This is a non-const routine that is public. + * + * The algorithm for this routine has undergone quite a bit of work. It probably needs more work. + * However, it seems now to be fairly robust. + * The key requirement is to find an initial pressure where both the liquid and the gas exist. This + * is not as easy as it sounds, and it gets exceedingly hard as the critical temperature is approached + * from below. + * Once we have this initial state, then we seek to equilibrate the gibbs free energies of the + * gas and liquid and use the formula + * + * dp = VdG + * + * to create an update condition for deltaP using + * + * - (Gliq - Ggas) = (Vliq - Vgas) (deltaP) + * + * + * + * @param TKelvin (input) Temperature (Kelvin) + * @param molarVolGas (return) Molar volume of the gas + * @param molarVolLiquid (return) Molar volume of the liquid + * + * @return Returns the saturation pressure at the given temperature + * + * @TODO Suggestions for the future would be to switch it to an algorithm that uses the gas molar volume + * and the liquid molar volumes as the fundamental unknowns. + * + */ + doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal &molarVolGas, + doublereal &molarVolLiquid) { + // we need this because this is a non-const routine that is public + setTemperature(TKelvin); + double tcrit = critTemperature(); + double RhoLiquid, RhoGas; + double RhoLiquidGood, RhoGasGood; + double densSave = density(); + double tempSave = temperature(); + double pres; + doublereal mw = meanMolecularWeight(); + bool conv = false; + if (TKelvin < tcrit) { + + pres = psatEst(TKelvin); + // trial value = Psat from correlation + int i; + doublereal volLiquid = liquidVolEst(TKelvin, pres); + RhoLiquidGood = mw / volLiquid; + RhoGasGood = pres * mw / (GasConstant * TKelvin); + doublereal delGRT, liqGRT, gasGRT; + int stab; + doublereal presLast = pres; + + +#ifdef DDDD + double pVec[100]; + int n = 0; + for (int i = 0; i < 50; i++) { + pVec[n++] = 3.40E6 + 0.01E5 * i; + } + + for (int i = 0; i < 50; i++) { + stab = corr0(TKelvin, pVec[i], RhoLiquid, RhoGas, liqGRT, gasGRT); + printf ("p = %g, T = %g, stab = %d, Rl = %g Rg = %g, Gl = %g, Gg = %g\n", + pVec[i], TKelvin, stab, RhoLiquid, RhoGas,liqGRT, gasGRT); + } +#endif + + /* + * First part of the calculation involves finding a pressure at which the + * gas and the liquid state coexists. + */ + doublereal presLiquid; + doublereal presGas; + doublereal presBase = pres; + bool foundLiquid = false; + bool foundGas = false; + + doublereal densLiquid = densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood); + if (densLiquid > 0.0) { + foundLiquid = true; + presLiquid = pres; + RhoLiquidGood = densLiquid; + } + if (!foundLiquid) { + for (int i = 0; i < 50; i++) { + pres = 1.1 * pres; + densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood); + if (densLiquid > 0.0) { + foundLiquid = true; + presLiquid = pres; + RhoLiquidGood = densLiquid; + break; + } + } + } + + pres = presBase; + doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood); + if (densGas <= 0.0) { + foundGas = false; + } else { + foundGas = true; + presGas = pres; + RhoGasGood = densGas; + } + if (!foundGas) { + for (int i = 0; i < 50; i++) { + pres = 0.9 * pres; + densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood); + if (densGas > 0.0) { + foundGas = true; + presGas = pres; + RhoGasGood = densGas; + break; + } + } + } + + if (foundGas && foundLiquid) { + if (presGas == presLiquid) { + pres = presGas; + goto startIteration; + } + pres = 0.5 * (presLiquid + presGas); + bool goodLiq; + bool goodGas; + for (int i = 0; i < 50; i++) { + + doublereal densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood); + if (densLiquid <= 0.0) { + goodLiq = false; + } else { + goodLiq = true; + RhoLiquidGood = densLiquid; + presLiquid = pres; + } + doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood); + if (densGas <= 0.0) { + goodGas = false; + } else { + goodGas = true; + RhoGasGood = densGas; + presGas = pres; + } + if (goodGas && goodLiq) { + break; + } + if (!goodLiq && !goodGas) { + pres = 0.5 * (pres + presLiquid); + } + if (goodLiq || goodGas) { + pres = 0.5 * (presLiquid + presGas); + } + + } + } + if (!foundGas || !foundLiquid) { + printf("error coundn't find a starting pressure\n"); + return (0.0); + } + if (presGas != presLiquid) { + printf("error coundn't find a starting pressure\n"); + return (0.0); + } + + startIteration: + pres = presGas; + presLast = pres; + RhoGas = RhoGasGood; + RhoLiquid = RhoLiquidGood; + + + /* + * Now that we have found a good pressure we can proceed with the algorithm. + */ + + for (i = 0; i < 20; i++) { + + stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT); + if (stab == 0) { + presLast = pres; + delGRT = liqGRT - gasGRT; + doublereal delV = mw * (1.0/RhoLiquid - 1.0/RhoGas); + doublereal dp = - delGRT * GasConstant * TKelvin / delV; + + if (fabs(dp) > 0.1 * pres) { + if (dp > 0.0) { + dp = 0.1 * pres; + } else { + dp = -0.1 * pres; + } + } + pres += dp; + + } else if (stab == -1) { + delGRT = 1.0E6; + if (presLast > pres) { + pres = 0.5 * (presLast + pres); + } else { + // we are stuck here - try this + pres = 1.1 * pres; + } + } else if (stab == -2) { + if (presLast < pres) { + pres = 0.5 * (presLast + pres); + } else { + // we are stuck here - try this + pres = 0.9 * pres; + } + } + molarVolGas = mw / RhoGas; + molarVolLiquid = mw / RhoLiquid; + + + if (fabs(delGRT) < 1.0E-8) { + conv = true; + break; + } + } + + molarVolGas = mw / RhoGas; + molarVolLiquid = mw / RhoLiquid; + // Put the fluid in the desired end condition + setState_TR(tempSave, densSave); + + return pres; + + + } else { + pres = critPressure(); + setState_TP(TKelvin, pres); + RhoGas = density(); + molarVolGas = mw / RhoGas; + molarVolLiquid = molarVolGas; + setState_TR(tempSave, densSave); + } + return pres; + } + + //==================================================================================================================== + // Calculate the pressure given the temperature and the molar volume + doublereal MixtureFugacityTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const { + throw CanteraError("MixtureFugacityTP::pressureCalc", "unimplemented"); + return 0.0; + } + //==================================================================================================================== + // Calculate the pressure given the temperature and the molar volume + doublereal MixtureFugacityTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const { + throw CanteraError("MixtureFugacityTP::dpdVCalc", "unimplemented"); + return 0.0; + } + //==================================================================================================================== + + /* + * void _updateStandardStateThermo() (protected, virtual, const) + * + * If m_useTmpStandardStateStorage is true, + * This function must be called for every call to functions in this + * class that need standard state properties. + * Child classes may require that it be called even if m_useTmpStandardStateStorage + * is not true. + * It checks to see whether the temperature has changed and + * thus the ss thermodynamics functions for all of the species + * must be recalculated. + * + * This + */ + void MixtureFugacityTP::_updateReferenceStateThermo() const { + double Tnow = temperature(); + + // If the temperature has changed since the last time these + // properties were computed, recompute them. + if (m_Tlast_ref != Tnow) { + m_spthermo->update(Tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]); + m_Tlast_ref = Tnow; + + // update the species Gibbs functions + int k; + for (k = 0; k < m_kk; k++) { + m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + } + doublereal pref = refPressure(); + if (pref <= 0.0) { + throw CanteraError("MixtureFugacityTP::_updateReferenceStateThermo()", "neg ref pressure"); + } + m_logc0 = log(pref/(GasConstant * Tnow)); + } + } + //==================================================================================================================== + + +} + + diff --git a/Cantera/src/thermo/MixtureFugacityTP.h b/Cantera/src/thermo/MixtureFugacityTP.h new file mode 100644 index 000000000..77e88b061 --- /dev/null +++ b/Cantera/src/thermo/MixtureFugacityTP.h @@ -0,0 +1,968 @@ +/** + * @file MixtureFugacityTP.h + * Header file for a derived class of ThermoPhase that handles + * non-ideal mixtures based on the fugacity models (see \ref thermoprops and + * class \link Cantera::MixtureFugacityTP MixtureFugacityTP\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: 2010-08-04 10:06:14 -0600 (Wed, 04 Aug 2010) $ + * $Revision: 547 $ + */ + +#ifndef CT_MIXTUREFUGACITYTP_H +#define CT_MIXTUREFUGACITYTP_H + +#include "ThermoPhase.h" +#include "VPSSMgr.h" +#include "ResidEval.h" + +namespace Cantera { + + class XML_Node; + class PDSS; + + + //! Various states of the Fugacity object. In general there can be multiple liquid + //! objects for a single phase identified with each species. + +#define FLUID_UNSTABLE -4 +#define FLUID_UNDEFINED -3 +#define FLUID_SUPERCRIT -2 +#define FLUID_GAS -1 +#define FLUID_LIQUID_0 0 +#define FLUID_LIQUID_1 1 +#define FLUID_LIQUID_2 2 +#define FLUID_LIQUID_3 3 +#define FLUID_LIQUID_4 4 +#define FLUID_LIQUID_5 5 +#define FLUID_LIQUID_6 6 +#define FLUID_LIQUID_7 7 +#define FLUID_LIQUID_8 8 +#define FLUID_LIQUID_9 9 + + + + /** + * @ingroup thermoprops + * + * This is a filter class for ThermoPhase that implements some prepatory + * steps for efficiently handling mixture of gases that whose standard states + * are defined as ideal gases, but which describe also non-ideal solutions. + * In addition a multicomponent liquid phase below the critical temperature of the + * mixture is also allowed. The main subclass will be a mixture Redlich-kwong class. + * + * Several concepts are introduced. The first concept is there are temporary + * variables for holding the species standard state values + * of Cp, H, S, G, and V at the last temperature and pressure called. These functions are not recalculated + * if a new call is made using the previous temperature and pressure. + * + * The other concept is that the current state of the mixture is tracked. + * The state variable is either GAS, LIQUID, or SUPERCRIT fluid. Additionally, + * the variable LiquidContent is used and may vary between 0 and 1. + * + * To support the above functionality, pressure and temperature variables, + * m_Plast_ss and m_Tlast_ss, are kept which store the last pressure and temperature + * used in the evaluation of standard state properties. + * + * Typically, only one liquid phase is allowed to be formed within these classes. + * Additionally, there is an inherent contradiction between three phase models and + * the ThermoPhase class. The ThermoPhase class is really only meant to represent a + * single instanteation of a phase. The three phase models may be in equilibrium with + * multiple phases of the fluid in equilibrium with each other. This has yet to be resolved. + * + * This class is usually used for non-ideal gases. + * + * + * @nosubgrouping + */ + class MixtureFugacityTP : public ThermoPhase { + + public: + + /*! + * + * @name Constructors and Duplicators for %MixtureFugacityTP + * + */ + //! Constructor. + MixtureFugacityTP(); + + //! Copy Constructor. + /*! + * @param b Object to be copied + */ + MixtureFugacityTP(const MixtureFugacityTP &b); + + //! Assignment operator + /*! + * @param b Object to be copied + */ + MixtureFugacityTP& operator=(const MixtureFugacityTP &b); + + //! Destructor. + virtual ~MixtureFugacityTP(); + + + //! Duplication routine + /*! + * @return Returns a duplication + */ + virtual ThermoPhase *duplMyselfAsThermoPhase() const; + + //@} + + /** + * @name Utilities (MixtureFugacityTP) + */ + //@{ + /** + * Equation of state type flag. The base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Constants defined for this purpose are + * listed in mix_defs.h. + */ + virtual int eosType() const { return 0; } + + //! This method returns the convention used in specification + //! of the standard state, of which there are currently two, + //! temperature based, and variable pressure based. + /*! + * Currently, there are two standard state conventions: + * - Temperature-based activities + * cSS_CONVENTION_TEMPERATURE 0 + * - default + * + * - Variable Pressure and Temperature -based activities + * cSS_CONVENTION_VPSS 1 + */ + virtual int standardStateConvention() const; + + //! Set the solution branch to force the ThermoPhase to exist on one branch or another + /*! + * @param solnBranch Branch that the solution is restricted to. + * the value -1 means gas. The value -2 means unrestricted. + * Values of zero or greater refer to species dominated condensed phases. + */ + virtual void setForcedSolutionBranch(int solnBranch); + + //! Report the solution branch which the solution is restricted to + /*! + * @return Branch that the solution is restricted to. + * the value -1 means gas. The value -2 means unrestricted. + * Values of zero or greater refer to species dominated condensed phases. + */ + virtual int forcedSolutionBranch() const; + + //! Report the solution branch which the solution is actually on + /*! + * @return Branch that the solution is restricted to. + * the value -1 means gas. The value -2 means superfluid.. + * Values of zero or greater refer to species dominated condensed phases. + */ + virtual int reportSolnBranchActual() const; + + + + //! Get the array of log concentration-like derivatives of the + //! log activity coefficients + /*! + * This function is a virtual method. For ideal mixtures + * (unity activity coefficients), this can return zero. + * Implementations should take the derivative of the + * logarithm of the activity coefficient with respect to the + * logarithm of the concentration-like variable (i.e. moles) + * that represents the standard state. + * This quantity is to be used in conjunction with derivatives of + * that concentration-like variable when the derivative of the chemical + * potential is taken. + * + * units = dimensionless + * + * @param dlnActCoeffdlnN_diag Output vector of derivatives of the + * log Activity Coefficients. length = m_kk + */ + virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const { + err("getdlnActCoeffdlnN_diag"); + } + + + //@} + /// @name Partial Molar Properties of the Solution (MixtureFugacityTP) + //@{ + + + //! Get the array of non-dimensional species chemical potentials + //! These are partial molar Gibbs free energies. + /*! + * \f$ \mu_k / \hat R T \f$. + * Units: unitless + * + * We close the loop on this function, here, calling + * getChemPotentials() and then dividing by RT. No need for child + * classes to handle. + * + * @param mu Output vector of non-dimensional species chemical potentials + * Length: m_kk. + */ + void getChemPotentials_RT(doublereal* mu) const; + + //@} + + /*! + * @name Properties of the Standard State of the Species in the Solution + * (MixtureFugacityTP) + * + * Within MixtureFugacityTP, these properties are calculated via a common routine, + * _updateStandardStateThermo(), + * which must be overloaded in inherited objects. + * The values are cached within this object, and are not recalculated unless + * the temperature or pressure changes. + */ + //@{ + + //! Get the array of chemical potentials at unit activity. + /*! + * These are the standard state chemical potentials \f$ \mu^0_k(T,P) + * \f$. The values are evaluated at the current temperature and pressure. + * + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure + * of the solution. + * + * @param mu Output vector of standard state chemical potentials. + * length = m_kk. units are J / kmol. + */ + virtual void getStandardChemPotentials(doublereal* mu) const; + + //! Get the nondimensional Enthalpy functions for the species + //! at their standard states at the current T and P of the solution. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure + * of the solution. + * + * @param hrt Output vector of standard state enthalpies. + * length = m_kk. units are unitless. + */ + virtual void getEnthalpy_RT(doublereal* hrt) const; + + + //! Get the array of nondimensional Enthalpy functions for the standard state species + /*! + * at the current T and P of the solution. + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure + * 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 + //! at their standard states of solution at the current T and P of the solution. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure + * 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 Gibbs functions for the standard + //! state of the species at the current T and P. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure + * of the solution. + * + * @param gpure Output vector of standard state + * Gibbs free energies. length = m_kk. + * units are J/kmol. + * + * @todo This could be eliminated. It doesn't fit into the current + * naming convention. + */ + void getPureGibbs(doublereal* gpure) const; + + //! Returns the vector of nondimensional internal Energies of the standard state at the current temperature + //! and pressure of the solution for each species. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure + * of the solution. + * + * \f[ + * u^{ss}_k(T,P) = h^{ss}_k(T) - P * V^{ss}_k + * \f] + * + * @param urt Output vector of nondimensional standard state + * internal energies. length = m_kk. + */ + virtual void getIntEnergy_RT(doublereal *urt) const; + + + //! Get the nondimensional Heat Capacities at constant + //! pressure for the standard state of the species at the current T and P. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure of the solution. + * + * @param cpr Output vector containing the + * the nondimensional Heat Capacities at constant + * pressure for the standard state of the species. + * Length: m_kk. + */ + virtual void getCp_R(doublereal* cpr) const; + + + //! Get the molar volumes of each species in their standard + //! states at the current T and P of the solution. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure of the solution. + * + * units = m^3 / kmol + * + * @param vol Output vector of species volumes. length = m_kk. + * units = m^3 / kmol + */ + virtual void getStandardVolumes(doublereal *vol) const; + + + //! Set the temperature of the phase + /*! + * Currently this passes down to setState_TP(). It does not + * make sense to calculate the standard state without first + * setting T and P. + * + * @param temp Temperature (kelvin) + */ + virtual void setTemperature(const doublereal temp); + + + //! Set the internally storred pressure (Pa) at constant + //! temperature and composition + /*! + * Currently this passes down to setState_TP(). It does not + * make sense to calculate the standard state without first + * setting T and P. + * + * @param p input Pressure (Pa) + */ + virtual void setPressure(doublereal p); + + + +protected: + /** + * Calculate the density of the mixture using the partial + * molar volumes and mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the pure + * species molar volumes. We have additionally specified + * in this class that the pure species molar volumes are + * independent of temperature and pressure. + * + * NOTE: This is a non-virtual function, which is not a + * member of the ThermoPhase base class. + */ + virtual void calcDensity(); + + public: + //! Set the temperature and pressure at the same time + /*! + * Note this function triggers a reevalulation of the standard + * state quantities. + * + * @param T temperature (kelvin) + * @param pres pressure (pascal) + */ + virtual void setState_TP(doublereal T, doublereal pres); + + //! Set the internally storred temperature (K) and density (kg/m^3) + /*! + * @param t Temperature in kelvin + * @param rho Density (kg/m^3) + */ + virtual void setState_TR(doublereal T, doublereal rho); + + //! Set the temperature (K), pressure (Pa), and mole fractions. + /*! + * Note, the mole fractions are set first before the pressure is set. + * Setting the pressure may involve the solution of a nonlinear equation. + * + * @param t Temperature (K) + * @param p Pressure (Pa) + * @param x Vector of mole fractions. + * Length is equal to m_kk. + */ + virtual void setState_TPX(doublereal t, doublereal p, const doublereal* x); + + + //! Set the mass fractions to the specified values, and then + //! normalize them so that they sum to 1.0. + /*! + * @param y Array of unnormalized mass fraction values (input). + * Must have a length greater than or equal to the number of species. + */ + virtual void setMassFractions(const doublereal* const y); + + + //!Set the mass fractions to the specified values without normalizing. + /*! + * This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set of + * equations. + * + * @param y Input vector of mass fractions. + * Length is m_kk. + */ + virtual void setMassFractions_NoNorm(const doublereal* const y); + + + + //! Set the mole fractions to the specified values, and then + //! normalize them so that they sum to 1.0. + /*! + * @param x Array of unnormalized mole fraction values (input). + * Must have a length greater than or equal to the number of species. + */ + virtual void setMoleFractions(const doublereal* const x); + + + //! Set the mole fractions to the specified values without normalizing. + /*! + * This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set ofequations. + * + * @param x Input vector of mole fractions. + * Length is m_kk. + */ + virtual void setMoleFractions_NoNorm(const doublereal* const x); + + + //! Set the concentrations to the specified values within the phase. + /*! + * @param c The input vector to this routine is in dimensional + * units. For volumetric phases c[k] is the + * concentration of the kth species in kmol/m3. + * For surface phases, c[k] is the concentration + * in kmol/m2. The length of the vector is the number + * of species in the phase. + */ + virtual void setConcentrations(const doublereal* const c); + + protected: + void setMoleFractions_NoState(const doublereal* const x); + + + public: + //! Returns the current pressure of the phase + /*! + * The pressure is an independent variable in this phase. Its current value + * is storred in the object MixtureFugacityTP. + * + * @return return the pressure in pascals. + */ + doublereal pressure() const { + return m_Pcurrent; + } + + + + + + protected: + + //! Updates the reference state thermodynamic functions at the current T of the solution. + /*! + * + * If m_useTmpStandardStateStorage is true, + * this function must be called for every call to functions in this + * class. It checks to see whether the temperature or pressure has changed and + * thus the ss thermodynamics functions for all of the species + * must be recalculated. + * + * This function is responsible for updating the following internal members, + * when m_useTmpStandardStateStorage is true. + * + * - m_hss_RT; + * - m_cpss_R; + * - m_gss_RT; + * - m_sss_R; + * - m_Vss + * + * If m_useTmpStandardStateStorage is not true, this function may be + * required to be called by child classes to update internal member data. + * + */ + virtual void _updateReferenceStateThermo() const; + public: + + //@} + /// @name Thermodynamic Values for the Species Reference States (MixtureFugacityTP) + /*! + * There are also temporary + * variables for holding the species reference-state values of Cp, H, S, and V at the + * last temperature and reference pressure called. These functions are not recalculated + * if a new call is made using the previous temperature. + * All calculations are done within the routine _updateRefStateThermo(). + */ + //@{ + + + //! Returns the vector of nondimensional + //! enthalpies of the reference state at the current temperature + //! of the solution and the reference pressure for the species. + /*! + * @param hrt Output vector contains the nondimensional enthalpies + * of the reference state of the species + * length = m_kk, units = dimensionless. + */ + virtual void getEnthalpy_RT_ref(doublereal *hrt) const; + +#ifdef H298MODIFY_CAPABILITY + //! Modify the value of the 298 K Heat of Formation of the standard state of + //! one species in the phase (J kmol-1) + /*! + * The 298K heat of formation is defined as the enthalpy change to create the standard state + * of the species from its constituent elements in their standard states at 298 K and 1 bar. + * + * @param k Index of the species + * @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar. + * units = J/kmol. + */ + void modifyOneHf298SS(const int k, const doublereal Hf298New); +#endif + + //! Returns the vector of nondimensional + //! Gibbs free energies of the reference state at the current temperature + //! of the solution and the reference pressure for the species. + /*! + * + * @param grt Output vector contains the nondimensional Gibbs free energies + * of the reference state of the species + * length = m_kk, units = dimensionless. + */ + virtual void getGibbs_RT_ref(doublereal *grt) const; + + protected: + //! Returns the vector of nondimensional + //! Gibbs free energies of the reference state at the current temperature + //! of the solution and the reference pressure for the species. + /*! + * @return Output vector contains the nondimensional Gibbs free energies + * of the reference state of the species + * length = m_kk, units = dimensionless. + */ + const vector_fp & gibbs_RT_ref() const; + + public: + /*! + * Returns the vector of the + * gibbs function of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * units = J/kmol + * + * @param g Output vector contain the Gibbs free energies + * of the reference state of the species + * length = m_kk, units = J/kmol. + */ + virtual void getGibbs_ref(doublereal *g) const; + + /*! + * Returns the vector of nondimensional + * entropies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * + * @param er Output vector contain the nondimensional entropies + * of the species in their reference states + * length: m_kk, units: dimensionless. + */ + virtual void getEntropy_R_ref(doublereal *er) const; + + /*! + * Returns the vector of nondimensional + * constant pressure heat capacities of the reference state + * at the current temperature of the solution + * and reference pressure for the species. + * + * @param cprt Output vector contains the nondimensional heat capacities + * of the species in their reference states + * length: m_kk, units: dimensionless. + */ + virtual void getCp_R_ref(doublereal *cprt) const; + + //! Get the molar volumes of the species reference states at the current + //! T and reference pressure of the solution. + /*! + * units = m^3 / kmol + * + * @param vol Output vector containing the standard state volumes. + * Length: m_kk. + */ + virtual void getStandardVolumes_ref(doublereal *vol) const; + + protected: + + + + //@} + + + public: + + //! @name Initialization Methods - For Internal use (VPStandardState) + /*! + * The following methods are used in the process of constructing + * the phase and setting its parameters from a specification in an + * input file. They are not normally used in application programs. + * To see how they are used, see files importCTML.cpp and + * ThermoFactory.cpp. + */ + //@{ + + + //! Set the initial state of the phase to the conditions specified in the state XML element. + /*! + * + * This method sets the temperature, pressure, and mole fraction vector to a set default value. + * + * @param state AN XML_Node object corresponding to + * the "state" entry for this phase in the input file. + */ + virtual void setStateFromXML(const XML_Node& state); + + //! @internal Initialize the object + /*! + * 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 after calling installSpecies() + * for each species in the phase. It's called before calling + * initThermoXML() for the phase. Therefore, it's the correct + * place for initializing vectors which have lengths equal to the + * number of species. + * + * @see importCTML.cpp + */ + virtual void initThermo(); + + //! Initialize a ThermoPhase object, potentially reading activity + //! coefficient information from an XML database. + /*! + * This routine initializes the lengths in the current object and + * then calls the parent routine. + * 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(). + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + virtual void initThermoXML(XML_Node& phaseNode, std::string id); + + + private: + //! @internal Initialize the internal lengths in this object. + /*! + * Note this is not a virtual function. + */ + void initLengths(); + + protected: + // Special Functions for fugacity classes + + + //! Calculate the value of z + /*! + * \f[ + * z = \frac{P v}{ R T} + * \f] + * + * returns the value of z + */ + doublereal z() const; + + //! Calculate the deviation terms for the total entropy of the mixture from the + //! ideal gas mixture + /* + * Here we use the current state conditions + * + * @return Returns the change in entropy in units of J kmol-1 K-1. + */ + virtual doublereal sresid() const; + + //! Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture + /* + * Here we use the current state conditions + * + * @return Returns the change in entropy in units of J kmol-1. + */ + virtual doublereal hresid() const; + + + //! Estimate for the saturation pressure + /*! + * Note: this is only used as a starting guess for later routines that actually calculate an + * accurate value for the saturation pressure. + * + * @param TKelvin temperature in kelvin + * + * @return returns the estimated saturation pressure at the given temperature + */ + virtual doublereal psatEst(doublereal TKelvin) const; + + //! Estimate for the molar volume of the liquid + /*! + * Note: this is only used as a starting guess for later routines that actually calculate an + * accurate value for the liquid molar volume. + * This routine doesn't change the state of the system. + * + * @param TKelvin temperature in kelvin + * @param pres Pressure in Pa. This is used as an initial guess. If the routine + * needs to change the pressure to find a stable liquid state, the + * new pressure is returned in this variable. + * + * @return Returns the estimate of the liquid volume. If the liquid can't be found, this + * routine returns -1. + */ + virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const; + + protected: + //! Calculates the density given the temperature and the pressure and a guess at the density. + /*! + * Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this. + * This is protected because it is called during setState_TP() routines. Infinite loops would result + * if it were not protected. + * + * -> why is this not const? + * + * parameters: + * @param TKelvin Temperature in Kelvin + * @param pressure Pressure in Pascals (Newton/m**2) + * @param phaseRequested int representing the phase whose density we are requesting. If we put + * a gas or liquid phase here, we will attempt to find a volume in that + * part of the volume space, only, in this routine. A value of FLUID_UNDEFINED + * means that we will accept anything. + * + * @param rhoguess Guessed density of the fluid. A value of -1.0 indicates that there + * is no guessed density + * + * + * @return We return the density of the fluid at the requested phase. If we have not found any + * acceptable density we return a -1. If we have found an accectable density at a + * different phase, we return a -2. + */ + virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phaseRequested, + doublereal rhoguess); + + //! Utility routine in the calculation of the saturation pressure + /*! + * Private routine + * + * @param TKelvin temperature (kelvin) + * @param pres pressure (Pascal) + * @param densLiq Output density of liquid + * @param densGas output density of gas + * @param delGRT output delGRT + */ + int corr0(doublereal TKelvin, doublereal pre, doublereal &densLiq, + doublereal &densGas, doublereal &liqGRT, doublereal &gasGRT); + public: + //! Returns the Phase State flag for the current state of the object + /*! + * @param checkState If true, this function does a complete check to see where + * in paramters space we are + * + * There are three values: + * WATER_GAS below the critical temperature but below the critical density + * WATER_LIQUID below the critical temperature but above the critical density + * WATER_SUPERCRIT above the critical temperature + */ + int phaseState(bool checkState = false) const ; + + //! Return the value of the density at the liquid spinodal point (on the liquid side) + //! for the current temperature. + /*! + * @return returns the density with units of kg m-3 + */ + virtual doublereal densSpinodalLiquid() const; + + + //! Return the value of the density at the gas spinodal point (on the gas side) + //! for the current temperature. + /*! + * @return returns the density with units of kg m-3 + */ + virtual doublereal densSpinodalGas() const; + + + + + public: + //! Calculate the saturation pressure at the current mixture content for the given temperature + /*! + * @param TKelvin (input) Temperature (Kelvin) + * @param molarVolGas (return) Molar volume of the gas + * @param molarVolLiquid (return) Molar volume of the liquid + * + * @return Returns the saturation pressure at the given temperature + */ + doublereal calculatePsat(doublereal TKelvin, doublereal &molarVolGas, + doublereal &molarVolLiquid); + protected: + //! Calculate the pressure given the temperature and the molar volume + /*! + * Calculate the pressure given the temperature and the molar volume + * + * @param TKelvin temperature in kelvin + * @param molarVol molar volume ( m3/kmol) + * + * @return Returns the pressure. + */ + virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const; + + + //! Calculate the pressure and the pressure derivative given the temperature and the molar volume + /*! + * Temperature and mole number are held constant + * + * @param TKelvin temperature in kelvin + * @param molarVol molar volume ( m3/kmol) + * + * @param presCalc Returns the pressure. + * + * @return Returns the derivative of the pressure wrt the molar volume + */ + virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const; + + + + virtual void updateMixingExpressions(); + + + //@} + + + class spinodalFunc : public Cantera::ResidEval + { + public: + + spinodalFunc(MixtureFugacityTP *tp); + + virtual int evalSS(const doublereal t, const doublereal * const y, doublereal * const r); + + MixtureFugacityTP *m_tp; + }; + + + protected: + + //! Current value of the pressurees + /*! + * Because the pressure is now a calculation, we store the result of the calculation whenever + * it is recalculated. + * + * units = Pascals + */ + doublereal m_Pcurrent; + + + //! Storage for the current values of the mole fractions of the species + /*! + * This vector is kept up-to-date when some the setState functions are called. + * + * The State object is allowed to com + * + * Therefore, it may be considered to be an independent variable. + * + + */ + std::vector moleFractions_; + + //! Current state of the fluid + /*! + * There are three possible states of the fluid + * FLUID_GAS + * FLUID_LIQUID + * FLUID_SUPERCRIT + */ + int iState_; + + + //! Force the system to be on a particular side of the spinodal curve + int forcedState_; + + //! The last temperature at which the reference state thermodynamic properties were calculated at. + mutable doublereal m_Tlast_ref; + + //! Temporary storage for log of p/rt + mutable doublereal m_logc0; + + //! Temporary storage for dimensionless reference state enthalpies + mutable array_fp m_h0_RT; + + //! Temporary storage for dimensionless reference state heat capacities + mutable array_fp m_cp0_R; + + //! Temporary storage for dimensionless reference state gibbs energies + mutable array_fp m_g0_RT; + + //! Temporary storage for dimensionless reference state entropies + mutable array_fp m_s0_R; + + spinodalFunc *fdpdv_; + private: + + //! MixtureFugacityTP has its own err routine + /*! + * @param msg Error message string + */ + doublereal err(std::string msg) const; + + }; +} + +#endif diff --git a/Cantera/src/thermo/NasaThermo.h b/Cantera/src/thermo/NasaThermo.h index 7d4d81fbd..95bbe4777 100644 --- a/Cantera/src/thermo/NasaThermo.h +++ b/Cantera/src/thermo/NasaThermo.h @@ -209,11 +209,12 @@ namespace Cantera { if (m_p0 < 0.0) { m_p0 = refPressure; } else if (fabs(m_p0 - refPressure) > 0.1) { - std::string logmsg = " WARNING NasaThermo: New Species, " + name + ", has a different reference pressure, " + std::string logmsg = " ERROR NasaThermo: New Species, " + name + ", has a different reference pressure, " + fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n"; writelog(logmsg); - logmsg = " This may become a fatal error in the future \n"; + logmsg = " This is now a fatal error\n"; writelog(logmsg); + throw CanteraError("install()", "species have different reference pressures"); } m_p0 = refPressure; } diff --git a/Cantera/src/thermo/RedlichKwongMFTP.cpp b/Cantera/src/thermo/RedlichKwongMFTP.cpp new file mode 100644 index 000000000..ec1627ba7 --- /dev/null +++ b/Cantera/src/thermo/RedlichKwongMFTP.cpp @@ -0,0 +1,1895 @@ +/** + * @file RedlichKwongMFTP.cpp + * Definition file for a derived class of ThermoPhase that assumes either + * an ideal gas or ideal solution approximation and handles + * variable pressure standard state methods for calculating + * thermodynamic properties (see \ref thermoprops and + * class \link Cantera::RedlichKwongMFTP RedlichKwongMFTP\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. + */ +/* + * $Author: hkmoffa $ + * $Date: 2009-11-09 16:36:49 -0700 (Mon, 09 Nov 2009) $ + * $Revision: 255 $ + */ + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "RedlichKwongMFTP.h" + + +#include "mix_defs.h" +#include "ThermoFactory.h" +#include "RootFind.h" + +using namespace std; + +namespace Cantera { +#ifdef WITH_REAL_GASSES + //==================================================================================================================== + /* + * Default constructor + */ + RedlichKwongMFTP::RedlichKwongMFTP() : + MixtureFugacityTP(), + m_standardMixingRules(0), + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + a_vec_Curr_(0), + b_vec_Curr_(0), + a_coeff_vec(0,0), + m_pc_Species(0), + m_tc_Species(0), + m_vc_Species(0), + NSolns_(0), + m_pp(0), + m_tmpV(0), + m_partialMolarVolumes(0), + dpdV_(0.0), + dpdT_(0.0), + dpdni_(0) + { + Vroot_[0] = 0.0; + Vroot_[1] = 0.0; + Vroot_[2] = 0.0; + } + //==================================================================================================================== + RedlichKwongMFTP::RedlichKwongMFTP(std::string infile, std::string id) : + MixtureFugacityTP(), + m_standardMixingRules(0), + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + a_vec_Curr_(0), + b_vec_Curr_(0), + a_coeff_vec(0,0), + m_pc_Species(0), + m_tc_Species(0), + m_vc_Species(0), + NSolns_(0), + m_pp(0), + m_tmpV(0), + m_partialMolarVolumes(0), + dpdV_(0.0), + dpdT_(0.0), + dpdni_(0) + { + Vroot_[0] = 0.0; + Vroot_[1] = 0.0; + Vroot_[2] = 0.0; + XML_Node* root = get_XML_File(infile); + if (id == "-") id = ""; + XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root); + if (!xphase) { + throw CanteraError("newPhase", + "Couldn't find phase named \"" + id + "\" in file, " + infile); + } + importPhase(*xphase, this); + } + //==================================================================================================================== + RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, std::string id) : + MixtureFugacityTP(), + m_standardMixingRules(0), + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + a_vec_Curr_(0), + b_vec_Curr_(0), + a_coeff_vec(0,0), + m_pc_Species(0), + m_tc_Species(0), + m_vc_Species(0), + NSolns_(0), + m_pp(0), + m_tmpV(0), + m_partialMolarVolumes(0), + dpdV_(0.0), + dpdT_(0.0), + dpdni_(0) + { + Vroot_[0] = 0.0; + Vroot_[1] = 0.0; + Vroot_[2] = 0.0; + XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, &phaseRefRoot); + if (!xphase) { + throw CanteraError("RedlichKwongMFTP::RedlichKwongMFTP()","Couldn't find phase named \"" + id + "\" in XML node"); + } + importPhase(*xphase, this); + } + + //==================================================================================================================== + RedlichKwongMFTP::RedlichKwongMFTP(int testProb) : + MixtureFugacityTP(), + m_standardMixingRules(0), + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + a_vec_Curr_(0), + b_vec_Curr_(0), + a_coeff_vec(0,0), + m_pc_Species(0), + m_tc_Species(0), + m_vc_Species(0), + NSolns_(0), + m_pp(0), + m_tmpV(0), + m_partialMolarVolumes(0), + dpdV_(0.0), + dpdT_(0.0), + dpdni_(0) + { + std::string infile = "co2_redlichkwong.xml"; + std::string id; + if (testProb == 1) { + infile = "co2_redlichkwong.xml"; + id = "carbondioxide"; + } else { + throw CanteraError("", "test prob = 1 only"); + } + XML_Node* root = get_XML_File(infile); + if (id == "-") id = ""; + XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root); + if (!xphase) { + throw CanteraError("newPhase", "Couldn't find phase named \"" + id + "\" in file, " + infile); + } + importPhase(*xphase, this); + } + //==================================================================================================================== + /* + * Copy Constructor: + * + * Note this stuff will not work until the underlying phase + * has a working copy constructor. + * + * The copy constructor just calls the assignment operator + * to do the heavy lifting. + */ + RedlichKwongMFTP::RedlichKwongMFTP(const RedlichKwongMFTP &b) : + MixtureFugacityTP(), + m_standardMixingRules(0), + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + a_vec_Curr_(0), + b_vec_Curr_(0), + a_coeff_vec(0,0), + m_pc_Species(0), + m_tc_Species(0), + m_vc_Species(0), + NSolns_(0), + m_pp(0), + m_tmpV(0), + m_partialMolarVolumes(0), + dpdV_(0.0), + dpdT_(0.0), + dpdni_(0) + { + *this = b; + } + + //==================================================================================================================== + /* + * operator=() + * + * Note this stuff will not work until the underlying phase + * has a working assignment operator + */ + RedlichKwongMFTP& RedlichKwongMFTP:: + operator=(const RedlichKwongMFTP &b) { + if (&b != this) { + /* + * Mostly, this is a passthrough to the underlying + * assignment operator for the ThermoPhae parent object. + */ + MixtureFugacityTP::operator=(b); + /* + * However, we have to handle data that we own. + */ + m_standardMixingRules = b.m_standardMixingRules; + m_formTempParam = b.m_formTempParam; + m_b_current = b.m_b_current; + m_a_current = b.m_a_current; + a_vec_Curr_ = b.a_vec_Curr_; + b_vec_Curr_ = b.b_vec_Curr_; + a_coeff_vec = b.a_coeff_vec; + + m_pc_Species = b.m_pc_Species; + m_tc_Species = b.m_tc_Species; + m_vc_Species = b.m_vc_Species; + NSolns_ = b.NSolns_; + Vroot_[0] = b.Vroot_[0]; + Vroot_[1] = b.Vroot_[1]; + Vroot_[2] = b.Vroot_[2]; + m_pp = b.m_pp; + m_tmpV = b.m_tmpV; + m_partialMolarVolumes = b.m_partialMolarVolumes; + dpdV_ = b.dpdV_; + dpdT_ = b.dpdT_; + dpdni_ = b.dpdni_; + } + return *this; + } + //==================================================================================================================== + /* + * ~RedlichKwongMFTP(): (virtual) + * + */ + RedlichKwongMFTP::~RedlichKwongMFTP() { + } + //==================================================================================================================== + /* + * Duplication function. + * This calls the copy constructor for this object. + */ + ThermoPhase* RedlichKwongMFTP::duplMyselfAsThermoPhase() const { + RedlichKwongMFTP* vptp = new RedlichKwongMFTP(*this); + return (ThermoPhase *) vptp; + } + //==================================================================================================================== + int RedlichKwongMFTP::eosType() const { + return cRedlichKwongMFTP; + } + + //==================================================================================================================== + /* + * ------------Molar Thermodynamic Properties ------------------------- + */ + //==================================================================================================================== + // Molar enthalpy. Units: J/kmol. + doublereal RedlichKwongMFTP::enthalpy_mole() const { + _updateReferenceStateThermo(); + doublereal rt = _RT(); + doublereal h_ideal = rt * mean_X(DATA_PTR(m_h0_RT)); + doublereal h_nonideal = hresid(); + return (h_ideal + h_nonideal); + } + //==================================================================================================================== + // Molar internal energy. Units: J/kmol. + doublereal RedlichKwongMFTP::intEnergy_mole() const { + doublereal p0 = pressure(); + doublereal md = molarDensity(); + return (enthalpy_mole() - p0 / md); + } + //==================================================================================================================== + // Molar entropy. Units: J/kmol/K. + doublereal RedlichKwongMFTP::entropy_mole() const { + _updateReferenceStateThermo(); + doublereal sr_ideal = GasConstant * (mean_X(DATA_PTR(m_s0_R)) + - sum_xlogx() - std::log(pressure()/m_spthermo->refPressure())); + doublereal sr_nonideal = sresid(); + return (sr_ideal + sr_nonideal); + } + //==================================================================================================================== + // Molar Gibbs function. Units: J/kmol. + doublereal RedlichKwongMFTP::gibbs_mole() const { + return enthalpy_mole() - temperature() * entropy_mole(); + } + //==================================================================================================================== + /// Molar heat capacity at constant pressure. Units: J/kmol/K. + doublereal RedlichKwongMFTP::cp_mole() const { + _updateReferenceStateThermo(); + doublereal TKelvin = temperature(); + doublereal sqt = sqrt(TKelvin); + doublereal mv = molarVolume(); + doublereal vpb = mv + m_b_current; + pressureDerivatives(); + doublereal cpref = GasConstant * mean_X(DATA_PTR(m_cp0_R)); + doublereal dadt = da_dt(); + doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0; + doublereal dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac + +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt)); + double cp = dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_; + return cp; + } + //==================================================================================================================== + /// Molar heat capacity at constant volume. Units: J/kmol/K. + doublereal RedlichKwongMFTP::cv_mole() const { + throw CanteraError("", "unimplemented"); + return cp_mole() - GasConstant; + } + //==================================================================================================================== + // Return the thermodynamic pressure (Pa). + /* + * Since the mass density, temperature, and mass fractions are stored, + * this method uses these values to implement the + * mechanical equation of state \f$ P(T, \rho, Y_1, \dots, Y_K) \f$. + * + * \f[ + * P = \frac{RT}{v-b_{mix}} - \frac{a_{mix}}{T^{0.5} v \left( v + b_{mix} \right) } + * \f] + * + */ + doublereal RedlichKwongMFTP::pressure() const { + + +#ifdef DEBUG_MODE + _updateReferenceStateThermo(); + + // Get a copy of the private variables stored in the State object + double rho = density(); + doublereal T = temperature(); + doublereal mmw = meanMolecularWeight(); + double molarV = mmw / rho; + + double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current)); + + if (fabs(pp -m_Pcurrent) > 1.0E-5 * m_Pcurrent) { + throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe"); + } +#endif + + return m_Pcurrent; + } + //==================================================================================================================== + void RedlichKwongMFTP::calcDensity() { + /* + * Calculate the molarVolume of the solution (m**3 kmol-1) + */ + + const doublereal * const dtmp = moleFractdivMMW(); + getPartialMolarVolumes(DATA_PTR(m_tmpV)); + double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp); + /* + * Set the density in the parent State object directly, + * by calling the State::setDensity() function. + */ + double dens = 1.0/invDens; + State::setDensity(dens); + + } + + //==================================================================================================================== + void RedlichKwongMFTP::setTemperature(const doublereal temp) { + State::setTemperature(temp); + _updateReferenceStateThermo(); + updateAB(); + } + //==================================================================================================================== + void RedlichKwongMFTP::setMassFractions(const doublereal * const x) { + MixtureFugacityTP::setMassFractions(x); + updateAB(); + } + //==================================================================================================================== + void RedlichKwongMFTP::setMassFractions_NoNorm(const doublereal * const x) { + MixtureFugacityTP::setMassFractions_NoNorm(x); + updateAB(); + } + //==================================================================================================================== + void RedlichKwongMFTP::setMoleFractions(const doublereal * const x) { + MixtureFugacityTP::setMoleFractions(x); + updateAB(); + } + //==================================================================================================================== + void RedlichKwongMFTP::setMoleFractions_NoNorm(const doublereal * const x) { + MixtureFugacityTP::setMoleFractions(x); + updateAB(); + } + //==================================================================================================================== + void RedlichKwongMFTP::setConcentrations(const doublereal * const c) { + MixtureFugacityTP::setConcentrations(c); + updateAB(); + } + + //==================================================================================================================== + doublereal RedlichKwongMFTP::isothermalCompressibility() const { + + + throw CanteraError("RedlichKwongMFTP::isothermalCompressibility() ", + "not implemented"); + + return 0.0; + } + //==================================================================================================================== + void RedlichKwongMFTP::getActivityConcentrations(doublereal* c) const { + + int k; + getPartialMolarVolumes(DATA_PTR(m_partialMolarVolumes)); + + for (k = 0; k < m_kk; k++) { + c[k] = moleFraction(k) / m_partialMolarVolumes[k]; + } + } + //==================================================================================================================== + /* + * Returns the standard concentration \f$ C^0_k \f$, which is used to normalize + * the generalized concentration. + */ + doublereal RedlichKwongMFTP::standardConcentration(int k) const { + + getStandardVolumes(DATA_PTR(m_tmpV)); + + return 1.0 / m_tmpV[k]; + + + + } + //==================================================================================================================== + /* + * Returns the natural logarithm of the standard + * concentration of the kth species + */ + doublereal RedlichKwongMFTP::logStandardConc(int k) const { + double c = standardConcentration(k); + double lc = std::log(c); + return lc; + } + //==================================================================================================================== + /* + * + * getUnitsStandardConcentration() + * + * Returns the units of the standard and general concentrations + * Note they have the same units, as their divisor 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 + * + * For EOS types other than cIdealSolidSolnPhase1, the default + * kmol/m3 holds for standard concentration units. For + * cIdealSolidSolnPhase0 type, the standard concentrtion is + * unitless. + */ + void RedlichKwongMFTP::getUnitsStandardConc(double *uA, int, int sizeUA) const { + //int eos = eosType(); + + for (int i = 0; i < sizeUA; i++) { + if (i == 0) uA[0] = 1.0; + if (i == 1) uA[1] = -nDim(); + if (i == 2) uA[2] = 0.0; + if (i == 3) uA[3] = 0.0; + if (i == 4) uA[4] = 0.0; + if (i == 5) uA[5] = 0.0; + } + + } + + //==================================================================================================================== + //! Get the array of non-dimensional activity coefficients at + //! the current solution temperature, pressure, and solution concentration. + /*! + * For ideal gases, the activity coefficients are all equal to one. + * + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + void RedlichKwongMFTP::getActivityCoefficients(doublereal *ac) const { + doublereal TKelvin = temperature(); + doublereal rt = TKelvin * GasConstant; + doublereal mv = molarVolume(); + doublereal sqt = sqrt(TKelvin); + doublereal vpb = mv + m_b_current; + doublereal vmb = mv - m_b_current; + + for (int k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; + } + } + doublereal pres = pressure(); + + for (int k = 0; k < m_kk; k++) { + ac[k] = ( - rt * log(pres * mv / rt) + + rt * log(mv / vmb) + + rt * b_vec_Curr_[k] / vmb + - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv) + + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) + - m_a_current / (m_b_current * sqt) * ( b_vec_Curr_[k]/vpb) + ); + } + for (int k = 0; k < m_kk; k++) { + ac[k] = exp(ac[k]/rt); + } + } + //==================================================================================================================== + /* + * ---- Partial Molar Properties of the Solution ----------------- + */ + //==================================================================================================================== + /* + * Get the array of non-dimensional species chemical potentials + * These are partial molar Gibbs free energies. + * \f$ \mu_k / \hat R T \f$. + * Units: unitless + * + * We close the loop on this function, here, calling + * getChemPotentials() and then dividing by RT. + */ + void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const{ + getChemPotentials(muRT); + doublereal invRT = 1.0 / _RT(); + for (int k = 0; k < m_kk; k++) { + muRT[k] *= invRT; + } + } + //==================================================================================================================== + void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const { + getGibbs_ref(mu); + doublereal xx; + doublereal rt = temperature() * GasConstant; + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(SmallNumber, moleFraction(k)); + mu[k] += rt*(log(xx)); + } + + doublereal TKelvin = temperature(); + doublereal mv = molarVolume(); + doublereal sqt = sqrt(TKelvin); + doublereal vpb = mv + m_b_current; + doublereal vmb = mv - m_b_current; + + for (int k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; + } + } + doublereal pres = pressure(); + doublereal refP = refPressure(); + + for (int k = 0; k < m_kk; k++) { + mu[k] += (rt * log(pres/refP) - rt * log(pres * mv / rt) + + rt * log(mv / vmb) + + rt * b_vec_Curr_[k] / vmb + - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv) + + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) + - m_a_current / (m_b_current * sqt) * ( b_vec_Curr_[k]/vpb) + ); + } + } + //==================================================================================================================== + void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const { + /* + * First we get the reference state contributions + */ + getEnthalpy_RT_ref(hbar); + doublereal rt = GasConstant * temperature(); + scale(hbar, hbar+m_kk, hbar, rt); + + /* + * We calculate dpdni_ + */ + doublereal TKelvin = temperature(); + doublereal mv = molarVolume(); + doublereal sqt = sqrt(TKelvin); + + doublereal vpb = mv + m_b_current; + doublereal vmb = mv - m_b_current; + + for (int k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; + } + } + + + + for (int k = 0; k < m_kk; k++) { + dpdni_[k] = rt/vmb + rt * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb) + + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb); + } + doublereal dadt = da_dt(); + doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0; + + for (int k = 0; k < m_kk; k++) { + m_tmpV[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter]; + } + } + + pressureDerivatives(); + doublereal fac2 = mv + TKelvin * dpdT_ / dpdV_; + + for (int k = 0; k < m_kk; k++) { + double hE_v = (mv * dpdni_[k] - rt - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac + + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k] + + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac); + hbar[k] = hbar[k] + hE_v; + + + hbar[k] -= fac2 * dpdni_[k]; + } + + } + //==================================================================================================================== + void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const { + getEntropy_R_ref(sbar); + doublereal r = GasConstant; + scale(sbar, sbar+m_kk, sbar, r); + doublereal TKelvin = temperature(); + doublereal sqt = sqrt(TKelvin); + doublereal mv = molarVolume(); + doublereal refP = refPressure(); + + for (int k = 0; k < m_kk; k++) { + doublereal xx = fmaxx(SmallNumber, moleFraction(k)); + sbar[k] += r * ( - log(xx)); + } + + for (int k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; + } + } + + for (int k = 0; k < m_kk; k++) { + m_tmpV[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter); + } + } + + + doublereal dadt = da_dt(); + doublereal fac = dadt - m_a_current / (2.0 * TKelvin); + doublereal vmb = mv - m_b_current; + doublereal vpb = mv + m_b_current; + + + for (int k = 0; k < m_kk; k++) { + sbar[k] -=( GasConstant * log(GasConstant * TKelvin / (refP * mv)) + + GasConstant + + GasConstant * log(mv/vmb) + + GasConstant * b_vec_Curr_[k]/vmb + + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv) + - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv) + + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac + - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac + ) ; + } + + pressureDerivatives(); + getPartialMolarVolumes(DATA_PTR(m_partialMolarVolumes)); + for (int k = 0; k < m_kk; k++) { + sbar[k] -= -m_partialMolarVolumes[k] * dpdT_; + } + } + //==================================================================================================================== + void RedlichKwongMFTP::getPartialMolarIntEnergies(doublereal* ubar) const { + getIntEnergy_RT(ubar); + doublereal rt = GasConstant * temperature(); + scale(ubar, ubar+m_kk, ubar, rt); + } + //==================================================================================================================== + void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const { + getCp_R(cpbar); + doublereal r = GasConstant; + scale(cpbar, cpbar+m_kk, cpbar, r); + } + //==================================================================================================================== + void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const { + // getStandardVolumes(vbar); + + + for (int k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; + } + } + + for (int k = 0; k < m_kk; k++) { + m_tmpV[k] = 0.0; + for (int i = 0; i < m_kk; i++) { + int counter = k + m_kk*i; + m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter); + } + } + + doublereal TKelvin = temperature(); + doublereal sqt = sqrt(TKelvin); + doublereal mv = molarVolume(); + + doublereal rt = GasConstant * TKelvin; + + doublereal vmb = mv - m_b_current; + doublereal vpb = mv + m_b_current; + + for (int k = 0; k < m_kk; k++) { + + doublereal num = (rt + rt * m_b_current/ vmb + rt * b_vec_Curr_[k] / vmb + + rt * m_b_current * b_vec_Curr_[k] /(vmb * vmb) + - 2.0 * m_pp[k] / (sqt * vpb) + + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb) + ); + + doublereal denom = (m_Pcurrent + rt * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb) + ); + + vbar[k] = num / denom; + } + + } + //==================================================================================================================== + doublereal RedlichKwongMFTP::critTemperature() const { + double pc, tc, vc; + double a0 = 0.0; + double aT = 0.0; + for (int i = 0; i < m_kk; i++) { + for (int j = 0; j 500.0) { + tmp2 = tmp / 500.; + tmp2 *= tmp2; + m_pp[k] = m_p0 * exp(500.) * tmp2; + } else { + m_pp[k] = m_p0 * exp(tmp); + } + pres += m_pp[k]; + } + // set state + setState_PX(pres, &m_pp[0]); + } + //==================================================================================================================== + /* + * Initialize the internal lengths. + * (this is not a virtual function) + */ + void RedlichKwongMFTP::initLengths() { + + + a_vec_Curr_.resize(m_kk * m_kk, 0.0); + b_vec_Curr_.resize(m_kk, 0.0); + + a_coeff_vec.resize(2, m_kk * m_kk, 0.0); + + + m_pc_Species.resize(m_kk, 0.0); + m_tc_Species.resize(m_kk, 0.0); + m_vc_Species.resize(m_kk, 0.0); + + + m_pp.resize(m_kk, 0.0); + m_tmpV.resize(m_kk, 0.0); + m_partialMolarVolumes.resize(m_kk, 0.0); + dpdni_.resize(m_kk, 0.0); + } + //==================================================================================================================== + /* + * Import and initialize a ThermoPhase object + * + * param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + * + * This routine initializes the lengths in the current object and + * then calls the parent routine. + */ + void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, std::string id) { + RedlichKwongMFTP::initLengths(); + + /* + * Check the model parameter for the Redlich-Kwong equation of state + * two are allowed + * RedlichKwong mixture of species, each of which are RK fluids + * RedlichKwongMFTP mixture of species with cross term coefficients + */ + if (phaseNode.hasChild("thermo")) { + XML_Node& thermoNode = phaseNode.child("thermo"); + std::string model = thermoNode["model"]; + if (model == "RedlichKwong" ) { + m_standardMixingRules = 1; + } else if (model == "RedlichKwongMFTP") { + m_standardMixingRules = 0; + } else { + throw CanteraError("RedlichKwongMFTP::initThermoXML", + "Unknown thermo model : " + model); + } + + + /* + * Go get all of the coefficients and factors in the + * activityCoefficients XML block + */ + XML_Node *acNodePtr = 0; + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& acNode = thermoNode.child("activityCoefficients"); + acNodePtr = &acNode; + int nC = acNode.nChildren(); + + /* + * Loop through the children getting multiple instances of + * parameters + */ + for (int i = 0; i < nC; i++) { + XML_Node &xmlACChild = acNodePtr->child(i); + string stemp = xmlACChild.name(); + string nodeName = lowercase(stemp); + /* + * Process a binary salt field, or any of the other XML fields + * that make up the Pitzer Database. Entries will be ignored + * if any of the species in the entry isn't in the solution. + */ + if (nodeName == "purefluidparameters") { + readXMLPureFluid(xmlACChild); + } + } + if (m_standardMixingRules == 1) { + applyStandardMixingRules(); + } + /* + * Loop through the children getting multiple instances of + * parameters + */ + for (int i = 0; i < nC; i++) { + XML_Node &xmlACChild = acNodePtr->child(i); + string stemp = xmlACChild.name(); + string nodeName = lowercase(stemp); + /* + * Process a binary salt field, or any of the other XML fields + * that make up the Pitzer Database. Entries will be ignored + * if any of the species in the entry isn't in the solution. + */ + if (nodeName == "crossfluidparameters") { + readXMLCrossFluid(xmlACChild); + } + } + + } + } + + for (int i = 0; i < m_kk; i++) { + double a0coeff = a_coeff_vec(0, i*m_kk + i); + double aTcoeff = a_coeff_vec(1, i*m_kk + i); + double ai = a0coeff + aTcoeff * 500.; + double bi = b_vec_Curr_[i]; + calcCriticalConditions(ai, bi, a0coeff, aTcoeff, m_pc_Species[i], m_tc_Species[i], m_vc_Species[i]); + } + + MixtureFugacityTP::initThermoXML(phaseNode, id); + } + //==================================================================================================================== + + void RedlichKwongMFTP::readXMLPureFluid(XML_Node &PureFluidParam) { + vector_fp vParams; + string xname = PureFluidParam.name(); + if (xname != "pureFluidParameters") { + throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", + "Incorrect name for processing this routine: " + xname); + } + + /* + * Read the species + * Find the index of the species in the current phase. It's not an error to not find the species + */ + string iName = PureFluidParam.attrib("species"); + if (iName == "") { + throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "no species attribute"); + } + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + int counter = iSpecies + m_kk * iSpecies; + int nParamsExpected, nParamsFound; + int num = PureFluidParam.nChildren(); + for (int iChild = 0; iChild < num; iChild++) { + XML_Node &xmlChild = PureFluidParam.child(iChild); + string stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + + if (nodeName == "a_coeff") { + string iModel = lowercase(xmlChild.attrib("model")); + if (iModel == "constant") { + nParamsExpected = 1; + } else if (iModel == "linear_a") { + nParamsExpected = 2; + if (m_formTempParam == 0) { + m_formTempParam = 1; + } + } else { + throw CanteraError("", "unknown model"); + } + + ctml::getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff"); + nParamsFound = vParams.size(); + if (nParamsFound != nParamsExpected) { + throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for a_coeff" + iName + ")", + "wrong number of params found"); + } + + for (int i = 0; i < nParamsFound; i++) { + a_coeff_vec(i, counter) = vParams[i]; + } + } else if (nodeName == "b_coeff") { + ctml::getFloatArray(xmlChild, vParams, true, "m3/kmol", "b_coeff"); + nParamsFound = vParams.size(); + if (nParamsFound != 1) { + throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for b_coeff" + iName + ")", + "wrong number of params found"); + } + b_vec_Curr_[iSpecies] = vParams[0]; + } + } + } + //==================================================================================================================== + void RedlichKwongMFTP::applyStandardMixingRules() { + int nParam = 2; + for (int i = 0; i < m_kk; i++) { + int icounter = i + m_kk * i; + for (int j = 0; j < m_kk; j++) { + if (i != j) { + int counter = i + m_kk * j; + int jcounter = j + m_kk * j; + for (int n = 0; n < nParam; n++) { + a_coeff_vec(n, counter) = sqrt(a_coeff_vec(n, icounter) * a_coeff_vec(n, jcounter)); + } + } + } + } + } + //==================================================================================================================== + + void RedlichKwongMFTP::readXMLCrossFluid(XML_Node &CrossFluidParam) { + vector_fp vParams; + string xname = CrossFluidParam.name(); + if (xname != "crossFluidParameters") { + throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", + "Incorrect name for processing this routine: " + xname); + } + + /* + * Read the species + * Find the index of the species in the current phase. It's not an error to not find the species + */ + string iName = CrossFluidParam.attrib("species1"); + if (iName == "") { + throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species1 attribute"); + } + int iSpecies = speciesIndex(iName); + if (iSpecies < 0) { + return; + } + string jName = CrossFluidParam.attrib("species2"); + if (iName == "") { + throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species2 attribute"); + } + int jSpecies = speciesIndex(jName); + if (jSpecies < 0) { + return; + } + + int counter = iSpecies + m_kk * jSpecies; + int counter0 = jSpecies + m_kk * iSpecies; + int nParamsExpected, nParamsFound; + int num = CrossFluidParam.nChildren(); + for (int iChild = 0; iChild < num; iChild++) { + XML_Node &xmlChild = CrossFluidParam.child(iChild); + string stemp = xmlChild.name(); + string nodeName = lowercase(stemp); + + if (nodeName == "a_coeff") { + string iModel = lowercase(xmlChild.attrib("model")); + if (iModel == "constant") { + nParamsExpected = 1; + } else if (iModel == "linear_a") { + nParamsExpected = 2; + if (m_formTempParam == 0) { + m_formTempParam = 1; + } + } else { + throw CanteraError("", "unknown model"); + } + + ctml::getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff"); + nParamsFound = vParams.size(); + if (nParamsFound != nParamsExpected) { + throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid(for a_coeff" + iName + ")", + "wrong number of params found"); + } + + for (int i = 0; i < nParamsFound; i++) { + a_coeff_vec(i, counter) = vParams[i]; + a_coeff_vec(i, counter0) = vParams[i]; + } + } + } + } + //==================================================================================================================== + void RedlichKwongMFTP::setParametersFromXML(const XML_Node& thermoNode) { + MixtureFugacityTP::setParametersFromXML(thermoNode); + std::string model = thermoNode["model"]; + + + } + //==================================================================================================================== + // Calculate the deviation terms for the total entropy of the mixture from the + // ideal gas mixture + /* + * Here we use the current state conditions + * + * @return Returns the change in entropy in units of J kmol-1 K-1. + */ + doublereal RedlichKwongMFTP::sresid() const{ + // note this agrees with tpx + doublereal rho = density(); + doublereal mmw = meanMolecularWeight(); + doublereal molarV = mmw / rho; + double hh = m_b_current / molarV; + doublereal zz = z(); + doublereal dadt = da_dt(); + doublereal T = temperature(); + doublereal sqT = sqrt(T); + doublereal fac = dadt - m_a_current / (2.0 * T); + double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current); + double sp = GasConstant * sresid_mol_R; + return sp; + } + //==================================================================================================================== + // Calculate the deviation terms for the total enthalpy of the mixture from the + // ideal gas mixture + /* + * Here we use the current state conditions + * + * @return Returns the change in entropy in units of J kmol-1. + */ + doublereal RedlichKwongMFTP::hresid() const{ + // note this agrees with tpx + doublereal rho = density(); + doublereal mmw = meanMolecularWeight(); + doublereal molarV = mmw / rho; + double hh = m_b_current / molarV; + doublereal zz = z(); + doublereal dadt = da_dt(); + doublereal T = temperature(); + doublereal sqT = sqrt(T); + doublereal fac = T * dadt - 3.0 *m_a_current / (2.0); + double hresid_mol = GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current); + return hresid_mol; + } + //==================================================================================================================== + // Estimate for the molar volume of the liquid + /* + * Note: this is only used as a starting guess for later routines that actually calculate an + * accurate value for the liquid molar volume. + * This routine doesn't change the state of the system. + * + * @param TKelvin temperature in kelvin + * @param pres Pressure in Pa. This is used as an initial guess. If the routine + * needs to change the pressure to find a stable liquid state, the + * new pressure is returned in this variable. + * + * @return Returns the estimate of the liquid volume. If the liquid can't be found, this + * routine returns -1. + */ + doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal &presGuess) const + { + + double v = m_b_current * 1.1; + double atmp; + double btmp; + calculateAB(TKelvin, atmp, btmp); + + doublereal pres = presGuess; + double pp = psatEst(TKelvin); + if (pres < pp) { + pres = pp; + } + double Vroot[3]; + +#ifdef NNN + if (TKelvin == 308.) { + double pVec[100]; + int n = 0; + for (int i = 0; i < 100; i++) { + pVec[n++] = 6.8E6 + 2.0E5 * i; + } + + for (int i = 0; i < 100; i++) { + int nsol = NicholsSolve(TKelvin, pVec[i], atmp, btmp, Vroot); + printf ("nsol = %d, p = %g, T = %g, v[0] = %g, v[1] %g, v[2] = %g\n", nsol, pVec[i], TKelvin, Vroot[0], Vroot[1], Vroot[2]); + } + } +#endif + + bool foundLiq = false; + int m = 0; + do { + + int nsol = NicholsSolve(TKelvin, pres, atmp, btmp, Vroot); + + // printf("nsol = %d\n", nsol); + // printf("liquidVolEst start: T = %g , p = %g, a = %g, b = %g\n", TKelvin, pres, m_a_current, m_b_current); + + if (nsol == 1 || nsol == 2) { + double pc = critPressure(); + if (pres > pc) { + foundLiq = true; + } + pres *= 1.04; + + } else { + foundLiq = true; + } + } while ((m < 100) && (!foundLiq)); + +#ifdef DONTUSE + int i; + double c; + double vnew; + double deltav; + double sqt = sqrt(TKelvin); + for (i = 0; i < 200; i++) { + c = bCalc * bCalc + bCalc * GasConstant * TKelvin / pres - atmp / (pres * sqt); + vnew = (1.0/c)*(v*v*v - GasConstant * TKelvin *v*v/pp - atmp * bCalc / (pres * sqt)); + deltav = vnew - v; + if (deltav > v*0.2) { + deltav = v * 0.2; + } else if (deltav < - (v * 0.2)) { + deltav = - v * 0.2; + } + v += deltav; + if (fabs(deltav) < 1.0E-6 * v) { + break; + } + } + if (i > 30) { + printf("liquidVolEst problem solve: T = %g , p = %g, a = %g, b = %g\n", TKelvin, pres, atmp, bCalc); + printf(" v final = %g\n", v); + } + if (fabs(deltav) > 1.0E-5 * v) { + throw CanteraError("RedlichKwongMFTP::liquidVolEst(T = " + fp2str(TKelvin) + ", " + fp2str(pres) + ")", + "failed to converge"); + } +#else + if (foundLiq) { + v = Vroot[0]; + presGuess = pres; + } else { + v = -1.0; + } +#endif + //printf (" RedlichKwongMFTP::liquidVolEst %g %g converged in %d its\n", TKelvin, pres, i); + return v; + } + //==================================================================================================================== + // Calculates the density given the temperature and the pressure and a guess at the density. + /* + * Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this. + * This is protected because it is called during setState_TP() routines. Infinite loops would result + * if it were not protected. + * + * -> why is this not const? + * + * parameters: + * @param TKelvin Temperature in Kelvin + * @param pressure Pressure in Pascals (Newton/m**2) + * @param phaseReqested int representing the phase whose density we are requesting. If we put + * a gas or liquid phase here, we will attempt to find a volume in that + * part of the volume space, only, in this routine. A value of FLUID_UNDEFINED + * means that we will accept anything. + * + * @param rhoguess Guessed density of the fluid. A value of -1.0 indicates that there + * is no guessed density + * + * + * @return We return the density of the fluid at the requested phase. If we have not found any + * acceptable density we return a -1. If we have found an accectable density at a + * different phase, we return a -2. + */ + doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess) { + + //setTemperature(TKelvin); + double tcrit = critTemperature(); + doublereal mmw = meanMolecularWeight(); + double densBase = 0.0; + if (rhoguess == -1.0) { + if (phaseRequested != FLUID_GAS) { + if (TKelvin > tcrit) { + rhoguess = presPa * mmw / (GasConstant * TKelvin); + } else { + if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) { + rhoguess = presPa * mmw / (GasConstant * TKelvin); + } else if (phaseRequested >= FLUID_LIQUID_0) { + double lqvol = liquidVolEst(TKelvin, presPa); + rhoguess = mmw / lqvol; + } + } + } else { + /* + * Assume the Gas phase initial guess, if nothing is + * specified to the routine + */ + rhoguess = presPa * mmw / (GasConstant * TKelvin); + } + + } + + + doublereal volguess = mmw / rhoguess; + NSolns_ = NicholsSolve(TKelvin, presPa, m_a_current, m_b_current, Vroot_); + + doublereal molarVolLast = Vroot_[0]; + if (NSolns_ >= 2) { + if (phaseRequested >= FLUID_LIQUID_0) { + molarVolLast = Vroot_[0]; + } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) { + molarVolLast = Vroot_[2]; + } else { + if (volguess > Vroot_[1]) { + molarVolLast = Vroot_[2]; + } else { + molarVolLast = Vroot_[0]; + } + } + } else if (NSolns_ == 1) { + if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) { + molarVolLast = Vroot_[0]; + } else { + //molarVolLast = Vroot_[0]; + //printf("DensityCalc(): Possible problem encountered\n"); + return -2.0; + } + } else if (NSolns_ == -1) { + if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) { + molarVolLast = Vroot_[0]; + } else if (TKelvin > tcrit) { + molarVolLast = Vroot_[0]; + } else { + // molarVolLast = Vroot_[0]; + //printf("DensityCalc(): Possible problem encountered\n"); + return -2.0; + } + } else { + molarVolLast = Vroot_[0]; + //printf("DensityCalc(): Possible problem encountered\n"); + return -1.0; + } + densBase = mmw / molarVolLast; + return densBase; + } + //==================================================================================================================== + // Return the value of the density at the liquid spinodal point (on the liquid side) + // for the current temperature. + /* + * @return returns the density with units of kg m-3 + */ + doublereal RedlichKwongMFTP::densSpinodalLiquid() const { + if (NSolns_ != 3) { + double dens = critDensity(); + return dens; + } + double vmax = Vroot_[1]; + double vmin = Vroot_[0]; + RootFind rf(fdpdv_); + rf.setPrintLvl(10); + rf.setTol(1.0E-5, 1.0E-10); + rf.setFuncIsGenerallyDecreasing(true); + + double vbest = 0.5 * (Vroot_[0]+Vroot_[1]); + double funcNeeded = 0.0; + + int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest); + if (status != ROOTFIND_SUCCESS) { + throw CanteraError(" RedlichKwongMFTP::densSpinodalLiquid() ", "didn't converge"); + } + doublereal mmw = meanMolecularWeight(); + doublereal rho = mmw / vbest; + return rho; + } + //==================================================================================================================== + // Return the value of the density at the gas spinodal point (on the gas side) + // for the current temperature. + /* + * @return returns the density with units of kg m-3 + */ + doublereal RedlichKwongMFTP::densSpinodalGas() const { + if (NSolns_ != 3) { + double dens = critDensity(); + return dens; + } + double vmax = Vroot_[2]; + double vmin = Vroot_[1]; + RootFind rf(fdpdv_); + rf.setPrintLvl(10); + rf.setTol(1.0E-5, 1.0E-10); + rf.setFuncIsGenerallyIncreasing(true); + + double vbest = 0.5 * (Vroot_[1]+Vroot_[2]); + double funcNeeded = 0.0; + + int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest); + if (status != ROOTFIND_SUCCESS) { + throw CanteraError(" RedlichKwongMFTP::densSpinodalGas() ", "didn't converge"); + } + doublereal mmw = meanMolecularWeight(); + doublereal rho = mmw / vbest; + return rho; + } + //==================================================================================================================== + // Calculate the pressure given the temperature and the molar volume + /* + * Calculate the pressure given the temperature and the molar volume + * + * @param TKelvin temperature in kelvin + * @param molarVol molar volume ( m3/kmol) + * + * @return Returns the pressure. + */ + doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const { + doublereal sqt = sqrt(TKelvin); + double pres = GasConstant * TKelvin / (molarVol - m_b_current) + - m_a_current / (sqt * molarVol * (molarVol + m_b_current)); + return pres; + } + //==================================================================================================================== + // Calculate the pressure and the pressure derivative given the temperature and the molar volume + /* + * Temperature and mole number are held constant + * + * @param TKelvin temperature in kelvin + * @param molarVol molar volume ( m3/kmol) + * + * @param presCalc Returns the pressure. + * + * @return Returns the derivative of the pressure wrt the molar volume + */ + doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const { + doublereal sqt = sqrt(TKelvin); + presCalc = GasConstant * TKelvin / (molarVol - m_b_current) + - m_a_current / (sqt * molarVol * (molarVol + m_b_current)); + + doublereal vpb = molarVol + m_b_current; + doublereal vmb = molarVol - m_b_current; + doublereal dpdv = (- GasConstant * TKelvin / (vmb * vmb) + + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb)); + return dpdv; + } + //==================================================================================================================== + + void RedlichKwongMFTP::pressureDerivatives() const { + doublereal TKelvin = temperature(); + doublereal mv = molarVolume(); + doublereal pres; + + dpdV_ = dpdVCalc(TKelvin, mv, pres); + + doublereal sqt = sqrt(TKelvin); + doublereal vpb = mv + m_b_current; + doublereal vmb = mv - m_b_current; + doublereal dadt = da_dt(); + doublereal fac = dadt - m_a_current/(2.0 * TKelvin); + + dpdT_ = (GasConstant / (vmb) - fac / (sqt * mv * vpb)); + } + //==================================================================================================================== + void RedlichKwongMFTP::updateMixingExpressions() { + updateAB(); + } + //==================================================================================================================== + void RedlichKwongMFTP::updateAB() { + double temp = temperature(); + int counter; + if (m_formTempParam == 1) { + for (int i = 0; i < m_kk; i++) { + for (int j = 0; j < m_kk; j++) { + counter = i * m_kk + j; + a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp; + } + } + } + + m_b_current = 0.0; + m_a_current = 0.0; + for (int i = 0; i < m_kk; i++) { + m_b_current += moleFractions_[i] * b_vec_Curr_[i]; + for (int j = 0; j < m_kk; j++) { + m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j]; + } + } + } + //==================================================================================================================== + void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal &aCalc, doublereal &bCalc) const { + int counter; + bCalc = 0.0; + aCalc = 0.0; + if (m_formTempParam == 1) { + for (int i = 0; i < m_kk; i++) { + bCalc += moleFractions_[i] * b_vec_Curr_[i]; + for (int j = 0; j < m_kk; j++) { + counter = i * m_kk + j; + doublereal a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp; + aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j]; + } + } + } else { + for (int i = 0; i < m_kk; i++) { + bCalc += moleFractions_[i] * b_vec_Curr_[i]; + for (int j = 0; j < m_kk; j++) { + counter = i * m_kk + j; + doublereal a_vec_Curr = a_coeff_vec(0,counter); + aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j]; + } + } + } + } + //==================================================================================================================== + doublereal RedlichKwongMFTP::da_dt() const { + + doublereal dadT = 0.0; + if (m_formTempParam == 1) { + for (int i = 0; i < m_kk; i++) { + for (int j = 0; j < m_kk; j++) { + int counter = i * m_kk + j; + dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j]; + } + } + } + return dadT; + } + //==================================================================================================================== + void RedlichKwongMFTP::calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff, + doublereal &pc, doublereal &tc, doublereal &vc) const { + if (m_formTempParam != 0) { + a = a0_coeff; + } + if (b <= 0.0) { + tc = 1000000.; + pc = 1.0E13; + vc = omega_vc * GasConstant * tc / pc; + return; + } + if (a <= 0.0) { + tc = 0.0; + pc = 0.0; + vc = 2.0 * b; + return; + } + double tmp = a * omega_b / (b * omega_a * GasConstant); + double pp = 2./3.; + doublereal sqrttc, f, dfdt, deltatc; + + if (m_formTempParam == 0) { + + tc = pow(tmp, pp); + } else { + tc = pow(tmp, pp); + for (int j = 0; j < 10; j++) { + sqrttc = sqrt(tc); + f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT_coeff * tc - a0_coeff; + dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT_coeff; + deltatc = - f / dfdt; + tc += deltatc; + } + if (deltatc > 0.1) { + throw CanteraError("RedlichKwongMFTP::calcCriticalConditions", "didn't converge"); + } + } + + pc = omega_b * GasConstant * tc / b; + vc = omega_vc * GasConstant * tc / pc; + } + + //==================================================================================================================== + // Solve the cubic equation of state + /* + * The R-K equation of state may be solved via the following formula + * + * V**3 - V**2(RT/P) - V(RTb/P - a/(P T**.5) + b*b) - (a b / (P T**.5)) = 0 + * + + * Returns the number of solutions found. If it only finds the liquid branch solution, it will return a -1 or a -2 + * instead of 1 or 2. If it returns 0, then there is an error. + * + */ + int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b, + doublereal Vroot[3]) const { + Vroot[0] = 0.0; + Vroot[1] = 0.0; + Vroot[2] = 0.0; + int nTurningPoints; + bool lotsOfNumError = false; + doublereal Vturn[2]; + if (TKelvin <= 0.0) { + throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "neg temperature"); + } + /* + * Derive the coefficients of the cubic polynomial to solve. + */ + doublereal an = 1.0; + doublereal bn = - GasConstant * TKelvin / pres; + doublereal sqt = sqrt(TKelvin); + doublereal cn = - (GasConstant * TKelvin * b / pres - a/(pres * sqt) + b * b); + doublereal dn = - (a * b / (pres * sqt)); + + double tmp = a * omega_b / (b * omega_a * GasConstant); + double pp = 2./3.; + double tc = pow(tmp, pp); + double pc = omega_b * GasConstant * tc / b; + double vc = omega_vc * GasConstant * tc / pc; + // Derive the center of the cubic, x_N + doublereal xN = - bn /(3 * an); + + + // Derive the value of delta**2. This is a key quantity that determines the number of turning points + doublereal delta2 = (bn * bn - 3 * an * cn) / (9 * an * an); + doublereal delta = 0.0; + + // Calculate a couple of ratios + doublereal ratio1 = 3.0 * an * cn / (bn * bn); + doublereal ratio2 = pres * b / (GasConstant * TKelvin); + if (fabs(ratio1) < 1.0E-7) { + //printf("NicholsSolve(): Alternative solution (p = %g T = %g)\n", pres, TKelvin); + doublereal ratio3 = a / (GasConstant * sqt) * pres / (GasConstant * TKelvin); + if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) { + doublereal z = 1.0; + for (int i = 0; i < 10; i++) { + doublereal znew = z / (z - ratio2) - ratio3 / (z + ratio1); + doublereal deltaz = znew - z; + z = znew; + if (fabs(deltaz) < 1.0E-14) { + break; + } + } + doublereal v = z * GasConstant * TKelvin / pres; + Vroot[0] = v; + return 1; + } + } + + + int nSolnValues; + nTurningPoints = 2; + +#ifdef PRINTPV + double V[100]; + int n = 0; + for (int i = 0; i < 90; i++) { + V[n++] = 0.030 + 0.005 * i; + } + double p1, presCalc; + for (int i = 0; i < n; i++) { + p1 = dpdVCalc(TKelvin, V[i], presCalc); + printf(" %13.5g %13.5g %13.5g \n", V[i], presCalc , p1); + } +#endif + + double h2 = 4. * an * an * delta2 * delta2 * delta2; + if (delta2 == 0.0) { + nTurningPoints = 1; + Vturn[0] = xN; + Vturn[1] = xN; + } else if (delta2 < 0.0) { + nTurningPoints = 0; + Vturn[0] = xN; + Vturn[1] = xN; + } else { + delta = sqrt(delta2); + Vturn[0] = xN - delta; + Vturn[1] = xN + delta; +#ifdef PRINTPV + double presCalc; + double p1 = dpdVCalc(TKelvin, Vturn[0], presCalc); + + double p2 = dpdVCalc(TKelvin, Vturn[1], presCalc); + + printf("p1 = %g p2 = %g \n", p1, p2); + p1 = dpdVCalc(TKelvin, 0.9*Vturn[0], presCalc); + printf("0.9 p1 = %g \n", p1); +#endif + } + + doublereal h = 2.0 * an * delta * delta2; + + doublereal yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn; + + doublereal desc = yN * yN - h2; + + if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) { + if (desc != 0.0) { + // this is for getting to other cases + printf("NicholsSolve(): numerical issues\n"); + throw CanteraError("NicholsSolve()", "numerical issues"); + } + desc = 0.0; + } + + if (desc < 0.0) { + nSolnValues = 3; + } else if (desc == 0.0) { + nSolnValues = 2; + // We are here as p goes to zero. + // double hleft = 3.0 * an * cn / (bn * bn); + //double ynleft = 9.0 * an * cn / (2.0 * bn * bn) - 27.0 * an * an * dn / (2.0 * bn * bn * bn); + //printf("hleft = %g , ynleft = %g\n", -3. / 2. * hleft, -ynleft); + //double h2left = - 3 * hleft + 3 * hleft * hleft - hleft * hleft * hleft; + //double y2left = - 2.0 * ynleft + ynleft * ynleft; + //printf("h2left = %g , yn2left = %g\n", h2left, y2left); + + } else if (desc > 0.0) { + nSolnValues = 1; + } + + /* + * One real root -> have to determine whether gas or liquid is the root + */ + if (desc > 0.0) { + doublereal tmpD = sqrt(desc); + doublereal tmp1 = (- yN + tmpD) / (2.0 * an); + doublereal sgn1 = 1.0; + if (tmp1 < 0.0) { + sgn1 = -1.0; + tmp1 = -tmp1; + } + doublereal tmp2 = (- yN - tmpD) / (2.0 * an); + doublereal sgn2 = 1.0; + if (tmp2 < 0.0) { + sgn2 = -1.0; + tmp2 = -tmp2; + } + doublereal p1 = pow(tmp1, 1./3.); + doublereal p2 = pow(tmp2, 1./3.); + + doublereal alpha = xN + sgn1 * p1 + sgn2 * p2; + Vroot[0] = alpha; + Vroot[1] = 0.0; + Vroot[2] = 0.0; + + double tmp = an * Vroot[0] * Vroot[0] * Vroot[0] + bn * Vroot[0] * Vroot[0] + cn * Vroot[0] + dn; + if (fabs (tmp) > 1.0E-4) { + lotsOfNumError = true; + } + + } else if (desc < 0.0) { + doublereal tmp = - yN/h; + + doublereal val = acos(tmp); + doublereal theta = val / 3.0; + + doublereal oo = 2. * Cantera::Pi / 3.; + doublereal alpha = xN + 2. * delta * cos(theta); + + doublereal beta = xN + 2. * delta * cos(theta + oo); + + doublereal gamma = xN + 2. * delta * cos(theta + 2.0 * oo); + + + Vroot[0] = beta; + Vroot[1] = gamma; + Vroot[2] = alpha; + + for (int i = 0; i < 3; i++) { + double tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; + if (fabs (tmp) > 1.0E-4) { + lotsOfNumError = true; + for (int j = 0; j < 3; j++) { + if (j != i) { + if (fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) { + writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + + fp2str(pres) + "): WARNING roots have merged: " + + fp2str(Vroot[i]) + ", " + fp2str(Vroot[j])); + writelogendl(); + } + } + } + } + } + } else if (desc == 0.0) { + if (yN == 0.0 && h == 0.0) { + Vroot[0] = xN; + Vroot[1] = xN; + Vroot[2] = xN; + } else { + // need to figure out whether delta is pos or neg + if (yN > 0.0) { + double tmp = pow(yN/(2*an), 1./3.); + if (fabs(tmp - delta) > 1.0E-9) { + throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected"); + } + Vroot[1] = xN + delta; + Vroot[0] = xN - 2.0*delta; // liquid phase root + } else { + double tmp = pow(yN/(2*an), 1./3.); + if (fabs(tmp - delta) > 1.0E-9) { + throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected"); + } + delta = -delta; + Vroot[0] = xN + delta; + Vroot[1] = xN - 2.0*delta; // gas phase root + } + } + for (int i = 0; i < 2; i++) { + double tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; + if (fabs (tmp) > 1.0E-4) { + lotsOfNumError = true; + } + } + } + + /* + * Unfortunately, there is a heavy amount of roundoff error due to bad conditioning in this + */ + double res, dresdV;; + for (int i = 0; i < nSolnValues; i++) { + for (int n = 0; n < 20; n++) { + res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; + if (fabs(res) < 1.0E-14) { + break; + } + dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn; + double del = - res / dresdV; + + Vroot[i] += del; + if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) { + break; + } + double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; + if (fabs(res2) < fabs(res)) { + continue; + } else { + Vroot[i] -= del; + Vroot[i] += 0.1 * del; + } + } + if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i])) ) { + writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + + fp2str(pres) + "): WARNING root didn't converge V = " + fp2str(Vroot[i]) ); + writelogendl(); + } + } + + if (nSolnValues == 1) { + if (TKelvin > tc) { + if (Vroot[0] < vc) { + nSolnValues = -1; + } + } else { + if (Vroot[0] < xN) { + nSolnValues = -1; + } + } + + } else { + if (nSolnValues == 2) { + if (delta > 0.0) { + nSolnValues = -2; + } + } + } + // writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + fp2str(pres) + "): finished"); + // writelogendl(); + return nSolnValues; + } + + +#endif +} + + diff --git a/Cantera/src/thermo/RedlichKwongMFTP.h b/Cantera/src/thermo/RedlichKwongMFTP.h new file mode 100644 index 000000000..8728fad52 --- /dev/null +++ b/Cantera/src/thermo/RedlichKwongMFTP.h @@ -0,0 +1,843 @@ +/** + * @file RedlichKwongMFTP.h + * Definition file for a derived class of ThermoPhase that assumes either + * an ideal gas or ideal solution approximation and handles + * variable pressure standard state methods for calculating + * thermodynamic properties (see \ref thermoprops and + * class \link Cantera::RedlichKwongMFTP RedlichKwongMFTP\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. + */ +/* + * $Author: hkmoffa $ + * $Date: 2009-11-09 16:36:49 -0700 (Mon, 09 Nov 2009) $ + * $Revision: 255 $ + */ + +#ifndef CT_REDLICHKWONGMFTP_H +#define CT_REDLICHKWONGMFTP_H + +#include "MixtureFugacityTP.h" + +namespace Cantera { + + class XML_Node; + class PDSS; + + /*! + * @name CONSTANTS - Models for the Standard State of IdealSolnPhase's + */ + //@{ + +#ifdef WITH_REAL_GASSES + + /** + * @ingroup thermoprops + * + * This class can handle either an ideal solution or an ideal gas approximation + * of a phase. + * + * + * @nosubgrouping + */ + class RedlichKwongMFTP : public MixtureFugacityTP { + + public: + + /*! + * + * @name Constructors and Duplicators for %RedlichKwongMFTP + * + */ + //! Base constructor. + RedlichKwongMFTP(); + + //! Construct and initialize a RedlichKwongMFTP ThermoPhase object + //! directly from an asci input file + /*! + * Working constructors + * + * The two constructors below are the normal way the phase initializes itself. They are shells that call + * the routine initThermo(), with a reference to the + * XML database to get the info for the phase. + * + * @param inputFile Name of the input file containing the phase XML data + * to set up the object + * @param id ID of the phase in the input file. Defaults to the empty string. + */ + RedlichKwongMFTP(std::string infile, std::string id=""); + + //! Construct and initialize a RedlichKwongMFTP ThermoPhase object + //! directly from an XML database + /*! + * @param phaseRef XML phase node containing the description of the phase + * @param id id attribute containing the name of the phase. (default is the empty string) + */ + RedlichKwongMFTP(XML_Node& phaseRef, std::string id = ""); + + //! This is a special constructor, used to replicate test problems + //! during the initial verification of the object + /*! + * + * test problems: + * 1: Pure CO2 problem + * input file = CO2_RedlickKwongMFTP.xml + * + * @param testProb Hard -coded test problem to instantiate. + * Current valid values are 1. + */ + RedlichKwongMFTP(int testProb); + + //! Copy Constructor + /*! + * Copy constructor for the object. Constructed object will be a clone of this object, but will + * also own all of its data. This is a wrapper around the assignment operator + * + * @param right Object to be copied. + */ + RedlichKwongMFTP(const RedlichKwongMFTP &right); + + //! Asignment operator + /*! + * Assignment operator for the object. Constructed object will be a clone of this object, but will + * also own all of its data. + * + * @param right Object to be copied. + */ + RedlichKwongMFTP& operator=(const RedlichKwongMFTP &right); + + //! Destructor. + virtual ~RedlichKwongMFTP(); + + + //! Duplicator from the ThermoPhase parent class + /*! + * Given a pointer to a ThermoPhase object, this function will + * duplicate the ThermoPhase object and all underlying structures. + * This is basically a wrapper around the copy constructor. + * + * @return returns a pointer to a ThermoPhase + */ + virtual ThermoPhase *duplMyselfAsThermoPhase() const; + + //@} + + /** + * @name Utilities (RedlichKwongMFTP) + */ + //@{ + /** + * Equation of state type flag. The base class returns + * zero. Subclasses should define this to return a unique + * non-zero value. Constants defined for this purpose are + * listed in mix_defs.h. + */ + virtual int eosType() const; + + //@} + + /// Molar enthalpy. Units: J/kmol. + virtual doublereal enthalpy_mole() const; + + /// Molar internal energy. Units: J/kmol. + virtual doublereal intEnergy_mole() const; + + /// Molar entropy. Units: J/kmol/K. + virtual doublereal entropy_mole() const; + + /// Molar Gibbs function. Units: J/kmol. + virtual doublereal gibbs_mole() const; + + /// Molar heat capacity at constant pressure. Units: J/kmol/K. + virtual doublereal cp_mole() const; + + /// Molar heat capacity at constant volume. Units: J/kmol/K. + virtual doublereal cv_mole() const; + + /** + * @} + * @name Mechanical Properties + * @{ + */ + + //! Return the thermodynamic pressure (Pa). + /*! + * Since the mass density, temperature, and mass fractions are stored, + * this method uses these values to implement the + * mechanical equation of state \f$ P(T, \rho, Y_1, \dots, Y_K) \f$. + * + * \f[ + * P = \frac{RT}{v-b_{mix}} - \frac{a_{mix}}{T^{0.5} v \left( v + b_{mix} \right) } + * \f] + * + */ + virtual doublereal pressure() const; + + //! 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; + + protected: + /** + * Calculate the density of the mixture using the partial + * molar volumes and mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the + * species standard state molar volumes. + * The species molar volumes may be functions + * of temperature and pressure. + * + * NOTE: This is a non-virtual function, which is not a + * member of the ThermoPhase base class. + */ + virtual void calcDensity(); + + //! Set the temperature (K) + /*! + * Overwritten setTemperature(double) from State.h. This + * function sets the temperature, and makes sure that + * the value propagates to underlying objects + * + * @todo Make State::setTemperature a virtual function + * + * @param temp Temperature in kelvin + */ + virtual void setTemperature(const doublereal temp); + + //! Set the mass fractions to the specified values, and then + //! normalize them so that they sum to 1.0. + /*! + * @param y Array of unnormalized mass fraction values (input). + * Must have a length greater than or equal to the number of species. + */ + virtual void setMassFractions(const doublereal* const y); + + //!Set the mass fractions to the specified values without normalizing. + /*! + * This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set of + * equations. + * + * @param y Input vector of mass fractions. + * Length is m_kk. + */ + virtual void setMassFractions_NoNorm(const doublereal* const y); + + //! Set the mole fractions to the specified values, and then + //! normalize them so that they sum to 1.0. + /*! + * @param x Array of unnormalized mole fraction values (input). + * Must have a length greater than or equal to the number of species. + */ + virtual void setMoleFractions(const doublereal* const x); + + //! Set the mole fractions to the specified values without normalizing. + /*! + * This is useful when the normalization + * condition is being handled by some other means, for example + * by a constraint equation as part of a larger set ofequations. + * + * @param x Input vector of mole fractions. + * Length is m_kk. + */ + virtual void setMoleFractions_NoNorm(const doublereal* const x); + + + //! Set the concentrations to the specified values within the phase. + /*! + * @param c The input vector to this routine is in dimensional + * units. For volumetric phases c[k] is the + * concentration of the kth species in kmol/m3. + * For surface phases, c[k] is the concentration + * in kmol/m2. The length of the vector is the number + * of species in the phase. + */ + virtual void setConcentrations(const doublereal* const c); + + + public: + + //! 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. Note that they may + * or may not have units of concentration --- they might be + * partial pressures, mole fractions, or surface coverages, + * for example. + * + * @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; + + //! Returns the standard concentration \f$ C^0_k \f$, which is used to normalize + //! the generalized concentration. + /*! + * This is defined as the concentration by which the generalized + * concentration is normalized to produce the activity. + * In many cases, this quantity will be the same for all species in a phase. + * Since the activity for an ideal gas mixture is + * simply the mole fraction, for an ideal gas \f$ C^0_k = P/\hat R T \f$. + * + * @param k Optional parameter indicating the species. The default + * is to assume this refers to species 0. + * @return + * Returns the standard Concentration in units of m3 kmol-1. + */ + virtual doublereal standardConcentration(int k=0) const; + + //! Returns the 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; + + //! 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 the default 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(double *uA, int k = 0, int sizeUA = 6) const; + + //! Get the array of non-dimensional activity coefficients at + //! the current solution temperature, pressure, and solution concentration. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure + * of the solution. The activities are based on this standard state. + * + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + virtual void getActivityCoefficients(doublereal* ac) const; + + + /// @name Partial Molar Properties of the Solution (RedlichKwongMFTP) + //@{ + + //! Get the array of non-dimensional species chemical potentials. + //! These are partial molar Gibbs free energies. + /*! + * \f$ \mu_k / \hat R T \f$. + * Units: unitless + * + * We close the loop on this function, here, calling + * getChemPotentials() and then dividing by RT. No need for child + * classes to handle. + * + * @param mu Output vector of non-dimensional species chemical potentials + * Length: m_kk. + */ + void getChemPotentials_RT(doublereal* mu) const; + + //! Get the species chemical potentials. Units: J/kmol. + /*! + * This function returns a vector of chemical potentials of the + * species in solution at the current temperature, pressure + * and mole fraction of the solution. + * + * @param mu Output vector of species chemical + * potentials. Length: m_kk. Units: J/kmol + */ + virtual void getChemPotentials(doublereal* mu) const; + + //! Get the species partial molar enthalpies. Units: J/kmol. + /*! + * @param hbar Output vector of species partial molar enthalpies. + * Length: m_kk. units are J/kmol. + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + //! Get the species partial molar entropies. Units: J/kmol/K. + /*! + * @param sbar Output vector of species partial molar entropies. + * Length = m_kk. units are J/kmol/K. + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + //! Get the species partial molar enthalpies. Units: J/kmol. + /*! + * @param ubar Output vector of speciar partial molar internal energies. + * Length = m_kk. units are J/kmol. + */ + virtual void getPartialMolarIntEnergies(doublereal* ubar) const; + + //! Get the partial molar heat capacities Units: J/kmol/K + /*! + * @param cpbar Output vector of species partial molar heat capacities + * at constant pressure. + * Length = m_kk. units are J/kmol/K. + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + //! Get the species partial molar volumes. Units: m^3/kmol. + /*! + * @param vbar Output vector of speciar partial molar volumes. + * Length = m_kk. units are m^3/kmol. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + //@} + + /*! + * @name Properties of the Standard State of the Species in the Solution + * + * Properties of the standard states are delegated to the VPSSMgr object. + * The values are cached within this object, and are not recalculated unless + * the temperature or pressure changes. + */ + //@{ + + //@} + + /// @name Thermodynamic Values for the Species Reference States (RedlichKwongMFTP) + /*! + * Properties of the reference states are delegated to the VPSSMgr object. + * The values are cached within this object, and are not recalculated unless + * the temperature or pressure changes. + */ + //@{ + + //@} + + + + //--------------------------------------------------------- + /// @name Critical State Properties. + /// These methods are only implemented by some subclasses, and may + /// be moved out of ThermoPhase at a later date. + + //@{ + + /// Critical temperature (K). + virtual doublereal critTemperature() const; + + /// Critical pressure (Pa). + virtual doublereal critPressure() const; + + /// Critical density (kg/m3). + virtual doublereal critDensity() const; + //@} + + + + + public: + + //! @name Initialization Methods - For Internal use (VPStandardState) + /*! + * The following methods are used in the process of constructing + * the phase and setting its parameters from a specification in an + * input file. They are not normally used in application programs. + * To see how they are used, see files importCTML.cpp and + * ThermoFactory.cpp. + */ + //@{ + + + //! 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. + * + * @param thermoNode An XML_Node object corresponding to + * the "thermo" entry for this phase in the input file. + */ + virtual void setParametersFromXML(const XML_Node& thermoNode); + + //! @internal Initialize the object + /*! + * 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(); + + + //!This method is used by the ChemEquil equilibrium solver. + /*! + * It sets the state such that the chemical potentials satisfy + * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m} + * \left(\frac{\lambda_m} {\hat R T}\right) \f] where + * \f$ \lambda_m \f$ is the element potential of element m. The + * temperature is unchanged. Any phase (ideal or not) that + * implements this method can be equilibrated by ChemEquil. + * + * @param lambda_RT Input vector of dimensionless element potentials + * The length is equal to nElements(). + */ + void setToEquilState(const doublereal* lambda_RT); + + //! Initialize a ThermoPhase object, potentially reading activity + //! coefficient information from an XML database. + /*! + * + * This routine initializes the lengths in the current object and + * then calls the parent routine. + * 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(). + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + virtual void initThermoXML(XML_Node& phaseNode, std::string id); + + private: + //! Read the pure species RedlichKwong input parameters + /*! + * @param pureFluidParam XML_Node for the pure fluid parameters + */ + void readXMLPureFluid(XML_Node &PureFluidParam); + + + //! Apply mixing rules for a coefficients + void applyStandardMixingRules(); + + + //! Read the cross species RedlichKwong input parameters + /*! + * @param pureFluidParam XML_Node for the cross fluid parameters + */ + void readXMLCrossFluid(XML_Node &PureFluidParam); + + + + //============================================================================== + private: + //! @internal Initialize the internal lengths in this object. + /*! + * Note this is not a virtual function and only handles + * this object + */ + void initLengths(); + + //============================================================================== + // Special functions inherited from MixtureFugacityTP + + protected: + + //! Calculate the deviation terms for the total entropy of the mixture from the + //! ideal gas mixture + /*! + * Here we use the current state conditions + * + * @return Returns the change in entropy in units of J kmol-1 K-1. + */ + virtual doublereal sresid() const; + + // Calculate the deviation terms for the total enthalpy of the mixture from the + // ideal gas mixture + /* + * Here we use the current state conditions + * + * @return Returns the change in enthalpy in units of J kmol-1. + */ + virtual doublereal hresid() const; + + //! Estimate for the molar volume of the liquid + /*! + * Note: this is only used as a starting guess for later routines that actually calculate an + * accurate value for the liquid molar volume. + * This routine doesn't change the state of the system. + * + * @param TKelvin temperature in kelvin + * @param pres Pressure in Pa. This is used as an initial guess. If the routine + * needs to change the pressure to find a stable liquid state, the + * new pressure is returned in this variable. + * + * @return Returns the estimate of the liquid volume. + */ + virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const; + + protected: + //! Calculates the density given the temperature and the pressure and a guess at the density. + /*! + * Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this. + * This is protected because it is called during setState_TP() routines. Infinite loops would result + * if it were not protected. + * + * -> why is this not const? + * + * parameters: + * @param TKelvin Temperature in Kelvin + * @param pressure Pressure in Pascals (Newton/m**2) + * @param phase int representing the phase whose density we are requesting. If we put + * a gas or liquid phase here, we will attempt to find a volume in that + * part of the volume space, only, in this routine. A value of FLUID_UNDEFINED + * means that we will accept anything. + * + * @param rhoguess Guessed density of the fluid. A value of -1.0 indicates that there + * is no guessed density + * + * + * @return We return the density of the fluid at the requested phase. If we have not found any + * acceptable density we return a -1. If we have found an accectable density at a + * different phase, we return a -2. + */ + virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phase, doublereal rhoguess); + + public: + //! Return the value of the density at the liquid spinodal point (on the liquid side) + //! for the current temperature. + /*! + * @return returns the density with units of kg m-3 + */ + virtual doublereal densSpinodalLiquid() const; + + + //! Return the value of the density at the gas spinodal point (on the gas side) + //! for the current temperature. + /*! + * @return returns the density with units of kg m-3 + */ + virtual doublereal densSpinodalGas() const; + + + + //! Calculate the pressure given the temperature and the molar volume + /*! + * Calculate the pressure given the temperature and the molar volume + * + * @param TKelvin temperature in kelvin + * @param molarVol molar volume ( m3/kmol) + * + * @return Returns the pressure. + */ + virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const; + + + //! Calculate the pressure and the pressure derivative given the temperature and the molar volume + /*! + * Temperature and mole number are held constant + * + * @param TKelvin temperature in kelvin + * @param molarVol molar volume ( m3/kmol) + * + * @param presCalc Returns the pressure. + * + * @return Returns the derivative of the pressure wrt the molar volume + */ + virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const; + + + //! Calculate dpdV and dpdT at the current conditions + /*! + * These are storred internally. + */ + void pressureDerivatives() const; + + + virtual void updateMixingExpressions(); + + + //! Update the a and b parameters + /*! + * The a and the b parameters depend on the mole fraction and the temperature. + * This function updates the internal numbers based on the state of the object. + */ + void updateAB(); + + + //! Calculate the a and the b parameters given the temperature + /*! + * + * This function doesn't change the internal state of the object, so it is a const + * function. It does use the storred mole fractions in the object. + * + * @param temp Temperature (TKelvin) + * + * @param aCalc (output) Returns the a value + * @param bCalc (output) Returns the b value. + */ + void calculateAB(doublereal temp, doublereal &aCalc, doublereal &bCalc) const; + + + //========================================================================================= + // Special functions not inherited from MixtureFugacityTP + + doublereal da_dt() const; + + void calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff, + doublereal &pc, doublereal &tc, doublereal &vc) const; + + + + int NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b, + doublereal Vroot[3]) const; + + //@} + //============================================================================== + protected: + + //! boolean indicating whether standard mixing rules are applied + /*! + * - 1 = Yes, there are standard cross terms in the a coefficient matrices. + * - 0 = No, there are nonstaandard cross terms in the a coefficient matrices. + */ + int m_standardMixingRules; + + //! Form of the temperature parameterization + /*! + * - 0 = There is no temperature parameterization of a or b + * - 1 = The a_ij parameter is a linear function of the temperature + */ + int m_formTempParam; + + + //! Value of b in the equation of state + /*! + * m_b is a function of the temperature and the mole fraction. + */ + doublereal m_b_current; + + //! Value of a in the equation of state + /*! + * a_b is a function of the temperature and the mole fraction. + */ + doublereal m_a_current; + + + vector_fp a_vec_Curr_; + vector_fp b_vec_Curr_; + + Array2D a_coeff_vec; + + + vector_fp m_pc_Species; + vector_fp m_tc_Species; + vector_fp m_vc_Species; + + int NSolns_; + + doublereal Vroot_[3]; + + + + //! Temporary storage - length = m_kk. + mutable vector_fp m_pp; + + //! Temporary storage - length = m_kk. + mutable vector_fp m_tmpV; + + // mutable vector_fp m_tmpV2; + + // Partial molar volumes of the species + mutable vector_fp m_partialMolarVolumes; + + + + //! The derivative of the pressure wrt the volume + /*! + * Calcualted at the current conditions + * temperature and mole number kept constant + */ + mutable doublereal dpdV_; + + //! The derivative of the pressure wrt the temperature + /*! + * Calcualted at the current conditions + * Total volume and mole number kept constant + */ + mutable doublereal dpdT_; + + //! Vector of derivatives of pressure wrt mole number + /*! + * Calcualted at the current conditions + * Total volume, temperature and other mole number kept constant + */ + mutable vector_fp dpdni_; + + public: + //! Omega constant for a -> value of a in terms of critical properties + /*! + * this was calculated from a small nonlinear solve + */ + static const doublereal omega_a = 4.27480233540E-01; + + //! Omega constant for b + static const doublereal omega_b = 8.66403499650E-02; + + //! Omega constant for the critical molar volume + static const doublereal omega_vc = 3.33333333333333E-01; + + + }; +#endif +} + +#endif diff --git a/Cantera/src/thermo/ShomateThermo.h b/Cantera/src/thermo/ShomateThermo.h index 06ef0bf2a..e82bbf8a4 100644 --- a/Cantera/src/thermo/ShomateThermo.h +++ b/Cantera/src/thermo/ShomateThermo.h @@ -201,12 +201,13 @@ namespace Cantera { if (m_p0 < 0.0) { m_p0 = refPressure; } else if (fabs(m_p0 - refPressure) > 0.1) { - std::string logmsg = " WARNING ShomateThermo: New Species, " + name + std::string logmsg = " ERROR ShomateThermo: New Species, " + name + ", has a different reference pressure, " + fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n"; writelog(logmsg); - logmsg = " This may become a fatal error in the future \n"; + logmsg = " This is now a fatal error\n"; writelog(logmsg); + throw CanteraError("install()", "Species have different reference pressures"); } m_p0 = refPressure; diff --git a/Cantera/src/thermo/SimpleThermo.h b/Cantera/src/thermo/SimpleThermo.h index cced50e0c..83bcb99e1 100644 --- a/Cantera/src/thermo/SimpleThermo.h +++ b/Cantera/src/thermo/SimpleThermo.h @@ -182,8 +182,9 @@ namespace Cantera { ", has a different reference pressure, " + fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n"; writelog(logmsg); - logmsg = " This may become a fatal error in the future \n"; + logmsg = " This is now a fatal error\n"; writelog(logmsg); + throw CanteraError("install()", "Species have different reference pressures"); } m_p0 = refPressure; } diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index a7181554c..9704cae20 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -36,6 +36,10 @@ #include "PureFluidPhase.h" #endif +#ifdef WITH_REAL_GASSES +#include "RedlichKwongMFTP.h" +#endif + #include "ConstDensityThermo.h" #include "SurfPhase.h" #include "EdgePhase.h" @@ -212,6 +216,13 @@ namespace Cantera { th = new PureFluidPhase; break; #endif + +#ifdef WITH_REAL_GASSES + case cRedlichKwongMFTP: + th = new RedlichKwongMFTP; + break; +#endif + #ifdef WITH_ELECTROLYTES case cHMW: th = new HMWSoln; diff --git a/Cantera/src/thermo/ThermoPhase.cpp b/Cantera/src/thermo/ThermoPhase.cpp index b014eaf89..746203f18 100644 --- a/Cantera/src/thermo/ThermoPhase.cpp +++ b/Cantera/src/thermo/ThermoPhase.cpp @@ -1012,7 +1012,7 @@ namespace Cantera { } return m_speciesData; } - + //==================================================================================================================== /* * Set the thermodynamic state. */ @@ -1038,7 +1038,7 @@ namespace Cantera { setDensity(rho); } } - + //==================================================================================================================== /* * Called by function 'equilibrate' in ChemEquil.h to transfer * the element potentials to this object after every successful diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h index d1b15ebbf..758919b3e 100644 --- a/Cantera/src/thermo/ThermoPhase.h +++ b/Cantera/src/thermo/ThermoPhase.h @@ -1544,7 +1544,7 @@ namespace Cantera { * @param x Vector of mole fractions. * Length is equal to m_kk. */ - void setState_TPX(doublereal t, doublereal p, const doublereal* x); + virtual void setState_TPX(doublereal t, doublereal p, const doublereal* x); //! Set the temperature (K), pressure (Pa), and mole fractions. /*! diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index b81981de7..3a4dec238 100644 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -73,6 +73,10 @@ namespace Cantera { const int cIdealSolnGasVPSS = 500; const int cIdealSolnGasVPSS_iscv = 501; + //! Fugacity Models + const int cMixtureFugacityTP = 700; + const int cRedlichKwongMFTP = 701; + const int cMargulesVPSSTP = 301; const int cPhaseCombo_Interaction = 305; diff --git a/config.h.in b/config.h.in index 767d95361..c13ec878b 100755 --- a/config.h.in +++ b/config.h.in @@ -183,6 +183,9 @@ typedef int ftnlen; // Fortran hidden string length type // models for electrolyte solutions. #undef WITH_ELECTROLYTES +// Enable Real Gasses +#undef WITH_REAL_GASSES + #undef WITH_PRIME // Enable the VCS NonIdeal equilibrium solver. This is diff --git a/configure b/configure index c746849d0..ed3d01168 100755 --- a/configure +++ b/configure @@ -2585,6 +2585,16 @@ _ACEOF fi +if test "$WITH_REAL_GASSES" = "y"; then + cat >>confdefs.h <<\_ACEOF +#define WITH_REAL_GASSES 1 +_ACEOF + + hdrs=$hdrs' RedlichKwongMFTP.h' + objs=$objs' RedlichKwongMFTP.o' +fi + + if test "$WITH_LATTICE_SOLID" = "y"; then cat >>confdefs.h <<\_ACEOF #define WITH_LATTICE_SOLID 1 @@ -10058,7 +10068,7 @@ fi # Provide some information about the compiler. -echo "$as_me:10061:" \ +echo "$as_me:10071:" \ "checking for Fortran 77 compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` { (eval echo "$as_me:$LINENO: \"$ac_compiler --version &5\"") >&5 @@ -10265,7 +10275,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_verb" -(eval echo $as_me:10268: \"$ac_link\") >&5 +(eval echo $as_me:10278: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS @@ -10343,7 +10353,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_cv_prog_f77_v" -(eval echo $as_me:10346: \"$ac_link\") >&5 +(eval echo $as_me:10356: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS diff --git a/configure.in b/configure.in index 5afbb896b..c58010bc6 100755 --- a/configure.in +++ b/configure.in @@ -486,6 +486,13 @@ if test "$WITH_PURE_FLUIDS" = "y"; then fi AC_SUBST(COMPILE_PURE_FLUIDS) +if test "$WITH_REAL_GASSES" = "y"; then + AC_DEFINE(WITH_REAL_GASSES) + hdrs=$hdrs' RedlichKwongMFTP.h' + objs=$objs' RedlichKwongMFTP.o' +fi + + if test "$WITH_LATTICE_SOLID" = "y"; then AC_DEFINE(WITH_LATTICE_SOLID) hdrs=$hdrs' LatticeSolidPhase.h' diff --git a/docs/install_examples/linux.64_sierra_gcc444_python264_numpy b/docs/install_examples/linux.64_sierra_gcc444_python264_numpy index 30e9e8181..32760d503 100755 --- a/docs/install_examples/linux.64_sierra_gcc444_python264_numpy +++ b/docs/install_examples/linux.64_sierra_gcc444_python264_numpy @@ -34,6 +34,9 @@ export WITH_VCSNONIDEAL WITH_H298MODIFY_CAPABILITY='y' export WITH_H298MODIFY_CAPABILITY +WITH_REAL_GASSES='y' +export WITH_REAL_GASSES + BUILD_MATLAB_TOOLBOX="y" export BUILD_MATLAB_TOOLBOX diff --git a/preconfig b/preconfig index 22259a98c..ebe22275d 100755 --- a/preconfig +++ b/preconfig @@ -228,6 +228,10 @@ WITH_IDEAL_SOLUTIONS=${WITH_IDEAL_SOLUTIONS:="y"} # models for electrolyte solutions WITH_ELECTROLYTES=${WITH_ELECTROLYTES:="y"} +# Enable Real Gas Equations of State +# This supports the multicomponent Redlich-Kwong equation of state. +WITH_REAL_GASSES=${WITH_REAL_GASSES:="n"} + # Enable generating phase models from PrIMe models. For more # information about PrIME, see http://www.primekinetics.org # WARNING: Support for PrIMe is experimental! diff --git a/test_problems/min_python/negATest/negATest_blessed.out b/test_problems/min_python/negATest/negATest_blessed.out index d170960c3..e63158035 100644 --- a/test_problems/min_python/negATest/negATest_blessed.out +++ b/test_problems/min_python/negATest/negATest_blessed.out @@ -1,7 +1,7 @@ Number of species = 12 Number of reactions = 12 rop [ 0:O ] = 0.447058 -rop [ 1:O2 ] = -0.0021443 +rop [ 1:O2 ] = -0.00214428 rop [ 2:N ] = 0 rop [ 3:NO ] = -279.361 rop [ 4:NO2 ] = 0.00214319 @@ -17,7 +17,7 @@ fwd_rop[ 0] = 479.304 rev_rop[ 0] = 97.94 fwd_rop[ 1] = -128.201 rev_rop[ 1] = -26.1964 fwd_rop[ 2] = 0 rev_rop[ 2] = 0 fwd_rop[ 3] = -0 rev_rop[ 3] = -0 -fwd_rop[ 4] = 0 rev_rop[ 4] = 1.10334e-06 +fwd_rop[ 4] = 0 rev_rop[ 4] = 1.08891e-06 fwd_rop[ 5] = 0 rev_rop[ 5] = 0 fwd_rop[ 6] = 0 rev_rop[ 6] = 0 fwd_rop[ 7] = 0 rev_rop[ 7] = 0