diff --git a/Cantera/src/thermo/Constituents.h b/Cantera/src/thermo/Constituents.h index a38622906..22028b964 100644 --- a/Cantera/src/thermo/Constituents.h +++ b/Cantera/src/thermo/Constituents.h @@ -280,14 +280,13 @@ namespace Cantera { doublereal charge = 0.0, doublereal size = 1.0); - //! Index of species named 'name' + //! Returns the index of a species named 'name' within the ThermoPhase /*! - * The first species added - * will have index 0, and the last one index nSpecies() - 1. + * The first species added will have index 0, and the last one index nSpecies() - 1. * * @param name String name of the species - * @return - * Returns the index of the species. + * @return Returns the index of the species. If the name is not found + * the value of -1 is returned. */ int speciesIndex(std::string name) const; diff --git a/Cantera/src/thermo/Elements.cpp b/Cantera/src/thermo/Elements.cpp index 1cd563d8e..2c9e61bdc 100644 --- a/Cantera/src/thermo/Elements.cpp +++ b/Cantera/src/thermo/Elements.cpp @@ -328,7 +328,6 @@ namespace Cantera { } - doublereal Elements::entropyElement298(int m) const { AssertThrowMsg(m_entropy298[m] != ENTROPY298_UNKNOWN, "Elements::entropy298", diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp index fd7e460c3..d0b14221f 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp @@ -23,6 +23,7 @@ #include "GibbsExcessVPSSTP.h" +#include using namespace std; namespace Cantera { @@ -65,13 +66,14 @@ namespace Cantera { moleFractions_ = b.moleFractions_; lnActCoeff_Scaled_ = b.lnActCoeff_Scaled_; dlnActCoeffdT_Scaled_ = b.dlnActCoeffdT_Scaled_; - dlnActCoeffdlnC_Scaled_ = b.dlnActCoeffdlnC_Scaled_; + dlnActCoeffdlnX_Scaled_ = b.dlnActCoeffdlnX_Scaled_; + dlnActCoeffdlnN_Scaled_ = b.dlnActCoeffdlnN_Scaled_; m_pp = b.m_pp; return *this; } - /** + /* * * ~GibbsExcessVPSSTP(): (virtual) * @@ -156,7 +158,9 @@ namespace Cantera { } void GibbsExcessVPSSTP::calcDensity() { - double *vbar = &m_pp[0]; + doublereal* vbar = NULL; + vbar = new doublereal[m_kk]; + // double *vbar = &m_pp[0]; getPartialMolarVolumes(vbar); doublereal vtotal = 0.0; @@ -165,6 +169,7 @@ namespace Cantera { } doublereal dd = meanMolecularWeight() / vtotal; State::setDensity(dd); + delete [] vbar; } void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p) { @@ -309,7 +314,6 @@ namespace Cantera { void GibbsExcessVPSSTP::initThermo() { initLengths(); VPStandardStateTP::initThermo(); - } @@ -320,110 +324,11 @@ namespace Cantera { moleFractions_.resize(m_kk); lnActCoeff_Scaled_.resize(m_kk); dlnActCoeffdT_Scaled_.resize(m_kk); - dlnActCoeffdlnC_Scaled_.resize(m_kk); + dlnActCoeffdlnX_Scaled_.resize(m_kk); + dlnActCoeffdlnN_Scaled_.resize(m_kk); m_pp.resize(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 GibbsExcessVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) { - - - VPStandardStateTP::initThermoXML(phaseNode, id); - } - /** - * Format a summary of the mixture state for output. - */ - std::string GibbsExcessVPSSTP::report(bool show_thermo) const { - - - char p[800]; - string s = ""; - try { - if (name() != "") { - sprintf(p, " \n %s:\n", name().c_str()); - s += p; - } - sprintf(p, " \n temperature %12.6g K\n", temperature()); - s += p; - sprintf(p, " pressure %12.6g Pa\n", pressure()); - s += p; - sprintf(p, " density %12.6g kg/m^3\n", density()); - s += p; - sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight()); - s += p; - - doublereal phi = electricPotential(); - sprintf(p, " potential %12.6g V\n", phi); - s += p; - - int kk = nSpecies(); - array_fp x(kk); - array_fp molal(kk); - array_fp mu(kk); - array_fp muss(kk); - array_fp acMolal(kk); - array_fp actMolal(kk); - getMoleFractions(&x[0]); - - getChemPotentials(&mu[0]); - getStandardChemPotentials(&muss[0]); - getActivities(&actMolal[0]); - - - if (show_thermo) { - sprintf(p, " \n"); - s += p; - sprintf(p, " 1 kg 1 kmol\n"); - s += p; - sprintf(p, " ----------- ------------\n"); - s += p; - sprintf(p, " enthalpy %12.6g %12.4g J\n", - enthalpy_mass(), enthalpy_mole()); - s += p; - sprintf(p, " internal energy %12.6g %12.4g J\n", - intEnergy_mass(), intEnergy_mole()); - s += p; - sprintf(p, " entropy %12.6g %12.4g J/K\n", - entropy_mass(), entropy_mole()); - s += p; - sprintf(p, " Gibbs function %12.6g %12.4g J\n", - gibbs_mass(), gibbs_mole()); - s += p; - sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n", - cp_mass(), cp_mole()); - s += p; - try { - sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n", - cv_mass(), cv_mole()); - s += p; - } - catch(CanteraError) { - sprintf(p, " heat capacity c_v \n"); - s += p; - } - } - - } catch (CanteraError) { - ; - } - return s; - } - } diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h index 185c2be27..a100e5f45 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.h +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h @@ -71,13 +71,14 @@ namespace Cantera { * \f$k\f$. * * GibbsExcessVPSSTP contains an internal vector with the current mole - * fraction vector. That's one of its primary usages. + * fraction vector. That's one of its primary usages. In order to keep the mole fraction + * vector constant, all of the setState functions are redesigned at this layer. * *

SetState Strategy

* - * The gibbsExcessVPSSTP object does not have a setState strategy. - * It's strictly an interfacial layer that writes the current mole fractions to the - * State object. + * All setState functions that set the internal state of the ThermoPhase object are + * overloaded at this level, so that a current mole fraction vector is maintained within + * the object. * * */ @@ -172,6 +173,7 @@ namespace Cantera { virtual void setPressure(doublereal p); protected: + /** * Calculate the density of the mixture using the partial * molar volumes and mole fractions as input @@ -301,6 +303,23 @@ namespace Cantera { virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) const { err("getdlnActCoeffdT"); } + + //! Get the array of change in the log activity coefficients w.r.t. change in state (change temp, change mole fractions) + /*! + * This function is a virtual class, but it first appears in GibbsExcessVPSSTP + * class and derived classes from GibbsExcessVPSSTP. + * + * This function is a virtual method. For ideal mixtures + * (unity activity coefficients), this can gradX/X. + * + * @param dT Input of temperature change + * @param dX Input vector of changes in mole fraction. length = m_kk + * @param dlnActCoeff Output vector of derivatives of the + * log Activity Coefficients. length = m_kk + */ + virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const { + err("getdlnActCoeff"); + } //! Get the array of log concentration-like derivatives of the //! log activity coefficients @@ -317,11 +336,33 @@ namespace Cantera { * * units = dimensionless * - * @param dlnActCoeffdlnC Output vector of derivatives of the + * @param dlnActCoeffdlnN Output vector of derivatives of the * log Activity Coefficients. length = m_kk */ - virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const { - err("getdlnActCoeffdlnC"); + virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const { + err("getdlnActCoeffdlnN"); + } + + //! 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. number of moles in + * in a unit volume. ) 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 dlnActCoeffdlnX Output vector of derivatives of the + * log Activity Coefficients. length = m_kk + */ + virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const { + err("getdlnActCoeffdlnX"); } //@} @@ -492,43 +533,14 @@ namespace Cantera { * @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); - - //! returns a summary of the state of the phase as a string - /*! - * @param show_thermo If true, extra information is printed out - * about the thermodynamic state of the system. - */ - virtual std::string report(bool show_thermo = true) const; - private: - //! Initialize lengths of local variables after all species have //! been identified. void initLengths(); - - - private: //! Error function /*! * Print an error string and exit @@ -548,6 +560,13 @@ namespace Cantera { protected: //! Storage for the current values of the mole fractions of the species + /*! + * This vector is kept up-to-date when the setState functions are called. + * Therefore, it may be considered to be an independent variable. + * + * Note in order to do this, the setState functions are redefined to always + * keep this vector current. + */ mutable std::vector moleFractions_; //! Storage for the current values of the activity coefficients of the @@ -562,7 +581,12 @@ namespace Cantera { //! Storage for the current derivative values of the //! gradients with respect to logarithm of the mole fraction of the //! log of theactivity coefficients of the species - mutable std::vector dlnActCoeffdlnC_Scaled_; + mutable std::vector dlnActCoeffdlnN_Scaled_; + + //! Storage for the current derivative values of the + //! gradients with respect to logarithm of the mole fraction of the + //! log of theactivity coefficients of the species + mutable std::vector dlnActCoeffdlnX_Scaled_; //! Temporary storage space that is fair game mutable std::vector m_pp; diff --git a/Cantera/src/thermo/IdealGasPhase.h b/Cantera/src/thermo/IdealGasPhase.h index b83cec970..915e7bfa0 100644 --- a/Cantera/src/thermo/IdealGasPhase.h +++ b/Cantera/src/thermo/IdealGasPhase.h @@ -330,6 +330,7 @@ namespace Cantera { */ IdealGasPhase(const IdealGasPhase &right); + //! Asignment operator /*! * Assignment operator for the object. Constructed diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp index df90bac16..5beb7725f 100644 --- a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -28,6 +28,7 @@ #include "mix_defs.h" #include +#include using namespace std; @@ -197,7 +198,8 @@ namespace Cantera { muNeutralMolecule_ = b.muNeutralMolecule_; gammaNeutralMolecule_ = b.gammaNeutralMolecule_; dlnActCoeffdT_NeutralMolecule_ = b.dlnActCoeffdT_NeutralMolecule_; - dlnActCoeffdlnC_NeutralMolecule_ = b.dlnActCoeffdlnC_NeutralMolecule_; + dlnActCoeffdlnX_NeutralMolecule_ = b.dlnActCoeffdlnX_NeutralMolecule_; + dlnActCoeffdlnN_NeutralMolecule_ = b.dlnActCoeffdlnN_NeutralMolecule_; return *this; } @@ -332,6 +334,14 @@ namespace Cantera { void IonsFromNeutralVPSSTP::getActivityConcentrations(doublereal* c) const { getActivities(c); } + + void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs,vector_fp& charges, std::vector& neutMolIndex){ + coeffs = fm_neutralMolec_ions_; + charges = m_speciesCharge; + neutMolIndex = fm_invert_ionForNeutral; + //for ( int k = 0; k < fm_neutralMolec_ions_[k]; k++ ) + // coeffs.push_back(fm_neutralMolec_ions_[k]); + } // Return the standard concentration for the kth species /* @@ -598,22 +608,47 @@ namespace Cantera { * * units = dimensionless * - * @param dlnActCoeffdlnC Output vector of log(mole fraction) + * @param dlnActCoeffdlnX Output vector of log(mole fraction) * derivatives of the log Activity Coefficients. * length = m_kk */ - void IonsFromNeutralVPSSTP::getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const { + void IonsFromNeutralVPSSTP::getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const { s_update_lnActCoeff(); - s_update_dlnActCoeff_dlnC(); + s_update_dlnActCoeff_dlnX(); for (int k = 0; k < m_kk; k++) { - dlnActCoeffdlnC[k] = dlnActCoeffdlnC_Scaled_[k]; + dlnActCoeffdlnX[k] = dlnActCoeffdlnX_Scaled_[k]; + } + } + + //! 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 Output vector of log(mole fraction) + * derivatives of the log Activity Coefficients. + * length = m_kk + */ + void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const { + s_update_lnActCoeff(); + s_update_dlnActCoeff_dlnN(); + + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnN[k] = dlnActCoeffdlnN_Scaled_[k]; } } - - - // This is temporary. We will get rid of this void IonsFromNeutralVPSSTP::setTemperature(const doublereal temp) { double p = pressure(); @@ -644,7 +679,10 @@ namespace Cantera { /* * Calculate the partial molar volumes, and then the density of the fluid */ - calcDensity(); + + //calcDensity(); + double dd = neutralMoleculePhase_->density(); + State::setDensity(dd); } // Calculate ion mole fractions from neutral molecule @@ -717,7 +755,7 @@ namespace Cantera { for (k = 0; k < m_kk; k++) { sum += moleFractions_[k]; } - if (fabs(sum) > 1.0E-11) { + if (fabs(sum) > 1.0E-11) { throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", "molefracts don't sum to one: " + fp2str(sum)); } @@ -791,7 +829,130 @@ namespace Cantera { sum += NeutralMolecMoleFractions_[k]; } for (k = 0; k < numNeutralMoleculeSpecies_; k++) { - NeutralMolecMoleFractions_[k] /= sum; + NeutralMolecMoleFractions_[k] /= sum; + } + + break; + + case cIonSolnType_SINGLECATION: + + throw CanteraError("eosType", "Unknown type"); + + break; + + case cIonSolnType_MULTICATIONANION: + + throw CanteraError("eosType", "Unknown type"); + break; + + default: + + throw CanteraError("eosType", "Unknown type"); + break; + + } + } + +// Calculate neutral molecule mole fractions + /* + * This routine calculates the neutral molecule mole + * fraction given the vector of ion mole fractions, + * i.e., the mole fractions from this ThermoPhase. + * Note, this routine basically assumes that there + * is charge neutrality. If there isn't, then it wouldn't + * make much sense. + * + * for the case of cIonSolnType_SINGLEANION, some slough + * in the charge neutrality is allowed. The cation number + * is followed, while the difference in charge neutrality + * is dumped into the anion mole number to fix the imbalance. + */ + void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal * const dx, doublereal *dy) const { + int k, icat, jNeut; + doublereal sumCat; + doublereal sumAnion; + doublereal fmij; + vector_fp y; + y.resize(numNeutralMoleculeSpecies_,0.0); + doublereal sumy, sumdy; + + //check sum dx = 0 + + //! Zero the vector we are trying to find. + for (k = 0; k < numNeutralMoleculeSpecies_; k++) { + dy[k] = 0.0; + } + + + // bool fmSimple = true; + + switch (ionSolnType_) { + + case cIonSolnType_PASSTHROUGH: + + for (k = 0; k < m_kk; k++) { + dy[k] = dx[k]; + } + break; + + case cIonSolnType_SINGLEANION: + + sumCat = 0.0; + sumAnion = 0.0; + + for (k = 0; k < (int) cationList_.size(); k++) { + //! Get the id for the next cation + icat = cationList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + if (jNeut >= 0) { + fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; + AssertTrace(fmij != 0.0); + dy[jNeut] += dx[icat] / fmij; + y[jNeut] += moleFractions_[icat] / fmij; + } + } + + for (k = 0; k < numPassThroughSpecies_; k++) { + icat = passThroughList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk]; + dy[jNeut] += dx[icat] / fmij; + y[jNeut] += moleFractions_[icat] / fmij; + } +#ifdef DEBUG_MODE_NOT +//check dy sum to zero + for (k = 0; k < m_kk; k++) { + moleFractionsTmp_[k] = dx[k]; + } + for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + for (k = 0; k < m_kk; k++) { + fmij = fm_neutralMolec_ions_[k + jNeut * m_kk]; + moleFractionsTmp_[k] -= fmij * dy[jNeut]; + } + } + for (k = 0; k < m_kk; k++) { + if (fabs(moleFractionsTmp_[k]) > 1.0E-13) { + //! Check to see if we have in fact found the inverse. + if (anionList_[0] != k) { + throw CanteraError("", "neutral molecule calc error"); + } else { + //! For the single anion case, we will allow some slippage + if (fabs(moleFractionsTmp_[k]) > 1.0E-5) { + throw CanteraError("", "neutral molecule calc error - anion"); + } + } + } + } +#endif + // Normalize the Neutral Molecule mole fractions + sumy = 0.0; + sumdy = 0.0; + for (k = 0; k < numNeutralMoleculeSpecies_; k++) { + sumy += y[k]; + sumdy += dy[k]; + } + for (k = 0; k < numNeutralMoleculeSpecies_; k++) { + dy[k] = dy[k]/sumy - y[k]*sumdy/sumy/sumy; } break; @@ -815,6 +976,7 @@ namespace Cantera { } } + void IonsFromNeutralVPSSTP::setMassFractions(const doublereal* const y) { GibbsExcessVPSSTP::setMassFractions(y); calcNeutralMoleculeMoleFractions(); @@ -836,7 +998,7 @@ namespace Cantera { void IonsFromNeutralVPSSTP::setMoleFractions_NoNorm(const doublereal* const x) { GibbsExcessVPSSTP::setMoleFractions_NoNorm(x); calcNeutralMoleculeMoleFractions(); - neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_)); + neutralMoleculePhase_->setMoleFractions_NoNorm(DATA_PTR(NeutralMolecMoleFractions_)); } @@ -1030,7 +1192,8 @@ namespace Cantera { muNeutralMolecule_.resize(numNeutralMoleculeSpecies_); gammaNeutralMolecule_.resize(numNeutralMoleculeSpecies_); dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); - dlnActCoeffdlnC_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); + dlnActCoeffdlnX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); + dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); } static double factorOverlap(const std::vector& elnamesVN , @@ -1136,9 +1299,12 @@ namespace Cantera { std::vector elemVectorI(nElementsI); vector fm_tmp(m_kk); - for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { - fm_invert_ionForNeutral[jNeut] = -1; + for (int k = 0; k < m_kk; k++) { + fm_invert_ionForNeutral[k] = -1; } + /* for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + fm_invert_ionForNeutral[jNeut] = -1; + }*/ for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { for (int m = 0; m < nElementsN; m++) { elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m); @@ -1185,12 +1351,16 @@ namespace Cantera { } bool notTaken = true; for (int iNeut = 0; iNeut < jNeut; iNeut++) { - if (fm_invert_ionForNeutral[iNeut] == k) { + if (fm_invert_ionForNeutral[k] == iNeut) { notTaken = false; } } if (notTaken) { - fm_invert_ionForNeutral[jNeut] = k; + fm_invert_ionForNeutral[k] = jNeut; + } + else{ + throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", + "Simple formula matrix generation failed, one cation is shared between two salts"); } } fm_neutralMolec_ions_[k + jNeut * m_kk] += fac; @@ -1243,7 +1413,7 @@ namespace Cantera { icat = cationList_[k]; jNeut = fm_invert_ionForNeutral[icat]; fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; - lnActCoeff_Scaled_[icat] = fmij * log(gammaNeutralMolecule_[jNeut]); + lnActCoeff_Scaled_[icat] = log(gammaNeutralMolecule_[jNeut])/fmij; } // Do the anion list @@ -1273,6 +1443,75 @@ namespace Cantera { } + + // get the gradient in the activity coefficients + + void IonsFromNeutralVPSSTP::getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const { + int k, icat, jNeut; + doublereal fmij; + int numNeutMolSpec; + /* + * Get the activity coefficients of the neutral molecules + */ + GibbsExcessVPSSTP *geThermo = dynamic_cast(neutralMoleculePhase_); + if (!geThermo) { + for ( k = 0; k < m_kk; k++ ){ + dlnActCoeff[k] = dX[k]/moleFractions_[k]; + } + return; + } + + numNeutMolSpec = geThermo->nSpecies(); + vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec); + vector_fp dX_NeutralMolecule(numNeutMolSpec); + + + getNeutralMoleculeMoleGrads(DATA_PTR(dX),DATA_PTR(dX_NeutralMolecule)); + + // All mole fractions returned to normal + + geThermo->getdlnActCoeff(dT, DATA_PTR(dX_NeutralMolecule), DATA_PTR(dlnActCoeff_NeutralMolecule)); + + switch (ionSolnType_) { + case cIonSolnType_PASSTHROUGH: + break; + case cIonSolnType_SINGLEANION: + + // Do the cation list + for (k = 0; k < (int) cationList_.size(); k++) { + //! Get the id for the next cation + icat = cationList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; + dlnActCoeff[icat] = dlnActCoeff_NeutralMolecule[jNeut]/fmij; + } + + // Do the anion list + icat = anionList_[0]; + jNeut = fm_invert_ionForNeutral[icat]; + dlnActCoeff[icat]= 0.0; + + // Do the list of neutral molecules + for (k = 0; k < numPassThroughSpecies_; k++) { + icat = passThroughList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + dlnActCoeff[icat] = dlnActCoeff_NeutralMolecule[jNeut]; + } + break; + + case cIonSolnType_SINGLECATION: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + case cIonSolnType_MULTICATIONANION: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + default: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + } + + } + // Update the temperatture derivative of the ln activity coefficients /* * This function will be called to update the internally storred @@ -1303,7 +1542,7 @@ namespace Cantera { icat = cationList_[k]; jNeut = fm_invert_ionForNeutral[icat]; fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; - dlnActCoeffdT_Scaled_[icat] = fmij * dlnActCoeffdT_NeutralMolecule_[jNeut]; + dlnActCoeffdT_Scaled_[icat] = dlnActCoeffdT_NeutralMolecule_[jNeut]/fmij; } // Do the anion list @@ -1336,7 +1575,7 @@ namespace Cantera { * This function will be called to update the internally storred * temperature derivative of the natural logarithm of the activity coefficients */ - void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnC() const { + void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX() const { int k, icat, jNeut; doublereal fmij; /* @@ -1344,11 +1583,11 @@ namespace Cantera { */ GibbsExcessVPSSTP *geThermo = dynamic_cast(neutralMoleculePhase_); if (!geThermo) { - fvo_zero_dbl_1(dlnActCoeffdlnC_Scaled_, m_kk); + fvo_zero_dbl_1(dlnActCoeffdlnX_Scaled_, m_kk); return; } - geThermo->getdlnActCoeffdlnC(DATA_PTR(dlnActCoeffdlnC_NeutralMolecule_)); + geThermo->getdlnActCoeffdlnX(DATA_PTR(dlnActCoeffdlnX_NeutralMolecule_)); switch (ionSolnType_) { case cIonSolnType_PASSTHROUGH: @@ -1361,19 +1600,19 @@ namespace Cantera { icat = cationList_[k]; jNeut = fm_invert_ionForNeutral[icat]; fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; - dlnActCoeffdlnC_Scaled_[icat] = fmij * dlnActCoeffdlnC_NeutralMolecule_[jNeut]; + dlnActCoeffdlnX_Scaled_[icat] = dlnActCoeffdlnX_NeutralMolecule_[jNeut]/fmij; } // Do the anion list icat = anionList_[0]; jNeut = fm_invert_ionForNeutral[icat]; - dlnActCoeffdT_Scaled_[icat]= 0.0; + dlnActCoeffdlnX_Scaled_[icat]= 0.0; // Do the list of neutral molecules for (k = 0; k < numPassThroughSpecies_; k++) { icat = passThroughList_[k]; jNeut = fm_invert_ionForNeutral[icat]; - dlnActCoeffdlnC_Scaled_[icat] = dlnActCoeffdlnC_NeutralMolecule_[jNeut]; + dlnActCoeffdlnX_Scaled_[icat] = dlnActCoeffdlnX_NeutralMolecule_[jNeut]; } break; @@ -1390,83 +1629,62 @@ namespace Cantera { } - /** - * Format a summary of the mixture state for output. - */ - std::string IonsFromNeutralVPSSTP::report(bool show_thermo) const { - - - char p[800]; - string s = ""; - try { - if (name() != "") { - sprintf(p, " \n %s:\n", name().c_str()); - s += p; - } - sprintf(p, " \n temperature %12.6g K\n", temperature()); - s += p; - sprintf(p, " pressure %12.6g Pa\n", pressure()); - s += p; - sprintf(p, " density %12.6g kg/m^3\n", density()); - s += p; - sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight()); - s += p; - - doublereal phi = electricPotential(); - sprintf(p, " potential %12.6g V\n", phi); - s += p; - - int kk = nSpecies(); - array_fp x(kk); - array_fp molal(kk); - array_fp mu(kk); - array_fp muss(kk); - array_fp acMolal(kk); - array_fp actMolal(kk); - getMoleFractions(&x[0]); - - getChemPotentials(&mu[0]); - getStandardChemPotentials(&muss[0]); - getActivities(&actMolal[0]); - - - if (show_thermo) { - sprintf(p, " \n"); - s += p; - sprintf(p, " 1 kg 1 kmol\n"); - s += p; - sprintf(p, " ----------- ------------\n"); - s += p; - sprintf(p, " enthalpy %12.6g %12.4g J\n", - enthalpy_mass(), enthalpy_mole()); - s += p; - sprintf(p, " internal energy %12.6g %12.4g J\n", - intEnergy_mass(), intEnergy_mole()); - s += p; - sprintf(p, " entropy %12.6g %12.4g J/K\n", - entropy_mass(), entropy_mole()); - s += p; - sprintf(p, " Gibbs function %12.6g %12.4g J\n", - gibbs_mass(), gibbs_mole()); - s += p; - sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n", - cp_mass(), cp_mole()); - s += p; - try { - sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n", - cv_mass(), cv_mole()); - s += p; - } - catch(CanteraError) { - sprintf(p, " heat capacity c_v \n"); - s += p; - } - } - - } catch (CanteraError) { - ; + /* + * This function will be called to update the internally storred + * temperature derivative of the natural logarithm of the activity coefficients + */ + void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const { + int k, icat, jNeut; + doublereal fmij; + /* + * Get the activity coefficients of the neutral molecules + */ + GibbsExcessVPSSTP *geThermo = dynamic_cast(neutralMoleculePhase_); + if (!geThermo) { + fvo_zero_dbl_1(dlnActCoeffdlnN_Scaled_, m_kk); + return; } - return s; + + geThermo->getdlnActCoeffdlnN(DATA_PTR(dlnActCoeffdlnN_NeutralMolecule_)); + + switch (ionSolnType_) { + case cIonSolnType_PASSTHROUGH: + break; + case cIonSolnType_SINGLEANION: + + // Do the cation list + for (k = 0; k < (int) cationList_.size(); k++) { + //! Get the id for the next cation + icat = cationList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; + dlnActCoeffdlnN_Scaled_[icat] = dlnActCoeffdlnN_NeutralMolecule_[jNeut]/fmij; + } + + // Do the anion list + icat = anionList_[0]; + jNeut = fm_invert_ionForNeutral[icat]; + dlnActCoeffdlnN_Scaled_[icat]= 0.0; + + // Do the list of neutral molecules + for (k = 0; k < numPassThroughSpecies_; k++) { + icat = passThroughList_[k]; + jNeut = fm_invert_ionForNeutral[icat]; + dlnActCoeffdlnN_Scaled_[icat] = dlnActCoeffdlnN_NeutralMolecule_[jNeut]; + } + break; + + case cIonSolnType_SINGLECATION: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + case cIonSolnType_MULTICATIONANION: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + default: + throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type"); + break; + } + } diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.h b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h index 88b05e786..2aa23319b 100644 --- a/Cantera/src/thermo/IonsFromNeutralVPSSTP.h +++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h @@ -404,6 +404,21 @@ namespace Cantera { */ virtual void getPartialMolarEntropies(doublereal* sbar) const; + //! Get the array of change in the log activity coefficients w.r.t. change in state (change temp, change mole fractions) + /*! + * This function is a virtual class, but it first appears in GibbsExcessVPSSTP + * class and derived classes from GibbsExcessVPSSTP. + * + * This function is a virtual method. For ideal mixtures + * (unity activity coefficients), this can gradX/X. + * + * @param dT Input of temperature change + * @param dX Input vector of changes in mole fraction. length = m_kk + * @param dlnActCoeff Output vector of derivatives of the + * log Activity Coefficients. length = m_kk + */ + virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const; + //! Get the array of log concentration-like derivatives of the //! log activity coefficients /*! @@ -411,19 +426,66 @@ namespace Cantera { * (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. + * logarithm of the concentration-like variable (i.e. mole fraction) + * 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 dlnActCoeffdlnC Output vector of log(mole fraction) + * @param dlnActCoeffdlnX Output vector of log(mole fraction) * derivatives of the log Activity Coefficients. * length = m_kk */ - virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const; + virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) 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. number of 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 Output vector of log(mole fraction) + * derivatives of the log Activity Coefficients. + * length = m_kk + */ + virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const; + + //! Get the Salt Dissociation Coefficients + //! Returns the vector of dissociation coefficients and vector of charges + virtual void getDissociationCoeffs(vector_fp& coeffs, vector_fp& charges, std::vector& neutMolIndex); + + virtual void getNeutralMolecMoleFractions(vector_fp& fracs){fracs=NeutralMolecMoleFractions_;} + + //! Calculate neutral molecule mole fractions + /*! + * This routine calculates the neutral molecule mole + * fraction given the vector of ion mole fractions, + * i.e., the mole fractions from this ThermoPhase. + * Note, this routine basically assumes that there + * is charge neutrality. If there isn't, then it wouldn't + * make much sense. + * + * for the case of cIonSolnType_SINGLEANION, some slough + * in the charge neutrality is allowed. The cation number + * is followed, while the difference in charge neutrality + * is dumped into the anion mole number to fix the imbalance. + */ + virtual void getNeutralMoleculeMoleGrads(const doublereal * const x, doublereal *y) const; + + virtual void getCationList(std::vector& cation){cation=cationList_;} + virtual void getAnionList(std::vector& anion){anion=anionList_;} + virtual void getSpeciesNames(std::vector& names){names=m_speciesNames;} //@} @@ -650,14 +712,6 @@ namespace Cantera { */ void initThermoXML(XML_Node& phaseNode, std::string id); - - //! returns a summary of the state of the phase as a string - /*! - * @param show_thermo If true, extra information is printed out - * about the thermodynamic state of the system. - */ - virtual std::string report(bool show_thermo = true) const; - private: @@ -680,6 +734,14 @@ namespace Cantera { */ void s_update_dlnActCoeffdT() const; + //! Update the change in the ln activity coefficients + /*! + * This function will be called to update the internally storred + * change of the natural logarithm of the activity coefficients + * w.r.t a change in state (temp, mole fraction, etc) + */ + void s_update_dlnActCoeff() const; + //! Update the derivative of the log of the activity coefficients //! wrt log(mole fraction) /*! @@ -687,7 +749,16 @@ namespace Cantera { * derivative of the natural logarithm of the activity coefficients * wrt logarithm of the mole fractions. */ - void s_update_dlnActCoeff_dlnC() const; + void s_update_dlnActCoeff_dlnX() const; + + //! Update the derivative of the log of the activity coefficients + //! wrt log(number of moles) + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the number of moles of given species. + */ + void s_update_dlnActCoeff_dlnN() const; private: @@ -726,7 +797,7 @@ namespace Cantera { //! Formula Matrix for composition of neutral molecules //! in terms of the molecules in this ThermoPhase /*! - * fm_neutralMolec_ions[ i + jNeut * NumNeut ] + * fm_neutralMolec_ions[ i + jNeut * m_kk ] * * This is the number of ions of type i in the neutral * molecule jNeut. @@ -735,6 +806,11 @@ namespace Cantera { //! Mapping between ion species and neutral molecule for quick invert. /*! + * + * fm_invert_ionForNeutral returns vector of int. Each element represents + * an ionic species and stores the value of the corresponding neutral + * molecule + * * For the case of fm_invert_simple_ = true, we assume that there * is a quick way to invert the formula matrix so that we can * quickly calculate the neutral molecule mole fraction @@ -818,8 +894,10 @@ namespace Cantera { mutable std::vector muNeutralMolecule_; mutable std::vector gammaNeutralMolecule_; + mutable std::vector dlnActCoeff_NeutralMolecule_; mutable std::vector dlnActCoeffdT_NeutralMolecule_; - mutable std::vector dlnActCoeffdlnC_NeutralMolecule_; + mutable std::vector dlnActCoeffdlnX_NeutralMolecule_; + mutable std::vector dlnActCoeffdlnN_NeutralMolecule_; }; diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp index 2017d8357..0d1ab06dc 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.cpp +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -19,7 +19,7 @@ #include "MargulesVPSSTP.h" #include "ThermoFactory.h" - +#include using namespace std; @@ -99,6 +99,12 @@ namespace Cantera { 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_; @@ -154,11 +160,20 @@ namespace Cantera { 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; @@ -166,6 +181,7 @@ namespace Cantera { 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) { @@ -217,7 +233,7 @@ namespace Cantera { */ void MargulesVPSSTP::constructPhaseFile(std::string inputFile, std::string id) { - if (inputFile.size() == 0) { + if ((int) inputFile.size() == 0) { throw CanteraError("MargulesVPSSTP:constructPhaseFile", "input file is null"); } @@ -274,7 +290,7 @@ namespace Cantera { */ void MargulesVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id) { string stemp; - if (id.size() > 0) { + if ((int) id.size() > 0) { string idp = phaseNode.id(); if (idp != id) { throw CanteraError("MargulesVPSSTP::constructPhaseXML", @@ -350,7 +366,7 @@ namespace Cantera { * take the exp of the internally storred coefficients. */ for (int k = 0; k < m_kk; k++) { - ac[k] = exp(lnActCoeff_Scaled_[k]); + ac[k] = exp(lnActCoeff_Scaled_[k]); } } @@ -472,7 +488,56 @@ namespace Cantera { } } - + /* + * ------------ 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 MargulesVPSSTP::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); + //cout << "species name(0) = " << speciesName(0) << endl; + //cout << "iA = " << speciesName(m_pSpecies_A_ij[0]) << endl; + //cout << "iB = " << speciesName(m_pSpecies_B_ij[0]) << endl; + + 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 MargulesVPSSTP::err(std::string msg) const { throw CanteraError("MargulesVPSSTP","Base class method " @@ -583,14 +648,54 @@ namespace Cantera { } + // 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_A - X_B)) + * he = X_A X_B(B + C X_B) */ + void MargulesVPSSTP::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; + //lnActCoeff_Scaled_[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1; + } + } + } + + + /* + // Not Right??? void MargulesVPSSTP::s_update_lnActCoeff() const { + int iA, iB; double XA, XB, g0 , g1; double T = temperature(); @@ -612,15 +717,52 @@ namespace Cantera { lnActCoeff_Scaled_[iB] += XA * XA * g0 + XA * XB * g1 * (2 * XA); } } + */ // 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_A - X_B)) + * he = X_A X_B(B + C X_B) */ void MargulesVPSSTP::s_update_dlnActCoeff_dT() const { + int iA, iB, iK, delAK, delBK; + double XA, XB, XK, g0 , g1; + double T = temperature(); + double RTT = GasConstant*T*T; + + fvo_zero_dbl_1(dlnActCoeffdT_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] / RTT; + g1 = -m_HE_c_ij[i] / RTT; + + dlnActCoeffdT_Scaled_[iK] += (delAK*XB+XA*delBK-XA*XB)*(g0+g1*XB)+XA*XB*(delBK-XB)*g1; + } + } + } + + /* Not Right??? + void MargulesVPSSTP::s_update_dlnActCoeff_dT() const {} + int iA, iB; doublereal XA, XB, h0 , h1; doublereal T = temperature(); @@ -642,6 +784,7 @@ namespace Cantera { dlnActCoeffdT_Scaled_[iB] += -(XA * XA * h0 + XA * XB * h1 * (2 * XA))/RTT; } } + */ void MargulesVPSSTP::getdlnActCoeffdT(doublereal *dlnActCoeffdT) const { s_update_dlnActCoeff_dT(); @@ -650,22 +793,123 @@ namespace Cantera { } } + // calculate the change of the log of the activity coefficients wrt change in state: dT, dX + /* + * This function will be called to calculate gradient of the + * logarithm of the activity coefficients based on gradients in temperature and mole fraction. + * + * he = X_A X_B(B + C X_B) + */ + void MargulesVPSSTP::getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal* dlnActCoeff) 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]; + dlnActCoeff[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 = dX[iA]; + dXB = dX[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; + + dlnActCoeff[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB + dlnActCoeffdT_Scaled_[iK]*dT; + } + } + } + // Update the derivative of the log of the activity coefficients wrt ln(X) /* * 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_A - X_B)) + * he = X_A X_B(B + C X_B) */ - void MargulesVPSSTP::s_update_dlnActCoeff_dlnC() const { + void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() 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_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; + + dlnActCoeffdlnN_Scaled_[iK] += 2*(delBK-XB) * (g0*(delAK-XA) + g1*(2*(delAK-XA)*XB + XA*(delBK-XB))); + } + dlnActCoeffdlnN_Scaled_[iK] = XK*dlnActCoeffdlnN_Scaled_[iK]-XK; + } + } + + void MargulesVPSSTP::s_update_dlnActCoeff_dlnX() const { + int iA, iB; doublereal XA, XB, g0 , g1; doublereal T = temperature(); - fvo_zero_dbl_1(dlnActCoeffdlnC_Scaled_, m_kk); + fvo_zero_dbl_1(dlnActCoeffdlnX_Scaled_, 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_Scaled_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB); + dlnActCoeffdlnX_Scaled_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB); + } + + /* + // Wrong!!! for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; iB = m_pSpecies_B_ij[i]; @@ -676,17 +920,25 @@ namespace Cantera { 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; - dlnActCoeffdlnC_Scaled_[iA] += XA * ( ( - 2.0 + 2.0 * XA ) * g0 + dlnActCoeffdlnX_Scaled_[iA] += XA * ( ( - 2.0 + 2.0 * XA ) * g0 + ( - 4.0 + 10.0 * XA - 6.0 * XA*XA ) * g1 ) ; - dlnActCoeffdlnC_Scaled_[iB] += XB * ( ( - 2.0 + 2.0 * XB ) * g0 + dlnActCoeffdlnX_Scaled_[iB] += XB * ( ( - 2.0 + 2.0 * XB ) * g0 + ( 2.0 - 8.0 * XB + 6.0 * XB*XB ) * g1 ) ; } + */ } + - void MargulesVPSSTP::getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const { - s_update_dlnActCoeff_dlnC(); + void MargulesVPSSTP::getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const { + s_update_dlnActCoeff_dlnN(); for (int k = 0; k < m_kk; k++) { - dlnActCoeffdlnC[k] = dlnActCoeffdlnC_Scaled_[k]; + dlnActCoeffdlnN[k] = dlnActCoeffdlnN_Scaled_[k]; + } + } + void MargulesVPSSTP::getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const { + s_update_dlnActCoeff_dlnX(); + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnX[k] = dlnActCoeffdlnX_Scaled_[k]; } } @@ -699,6 +951,12 @@ namespace Cantera { 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); @@ -769,7 +1027,7 @@ namespace Cantera { /* * Get the string containing all of the values */ - getFloatArray(xmlChild, vParams, true, "", "excessEnthalpy"); + getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy"); nParamsFound = vParams.size(); if (nParamsFound != 2) { @@ -785,7 +1043,7 @@ namespace Cantera { /* * Get the string containing all of the values */ - getFloatArray(xmlChild, vParams, true, "", "excessEntropy"); + getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy"); nParamsFound = vParams.size(); if (nParamsFound != 2) { @@ -797,90 +1055,42 @@ namespace Cantera { m_SE_c_ij[iSpot] = vParams[1]; } + if (nodeName == "excessvolume_enthalpy") { + /* + * Get the string containing all of the values + */ + getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MargulesVPSSTP::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 + */ + getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy"); + nParamsFound = vParams.size(); + + if (nParamsFound != 2) { + throw CanteraError("MargulesVPSSTP::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]; + } + } } - - /** - * Format a summary of the mixture state for output. - */ - std::string MargulesVPSSTP::report(bool show_thermo) const { - - - char p[800]; - string s = ""; - try { - if (name() != "") { - sprintf(p, " \n %s:\n", name().c_str()); - s += p; - } - sprintf(p, " \n temperature %12.6g K\n", temperature()); - s += p; - sprintf(p, " pressure %12.6g Pa\n", pressure()); - s += p; - sprintf(p, " density %12.6g kg/m^3\n", density()); - s += p; - sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight()); - s += p; - - doublereal phi = electricPotential(); - sprintf(p, " potential %12.6g V\n", phi); - s += p; - - int kk = nSpecies(); - array_fp x(kk); - array_fp molal(kk); - array_fp mu(kk); - array_fp muss(kk); - array_fp acMolal(kk); - array_fp actMolal(kk); - getMoleFractions(&x[0]); - - getChemPotentials(&mu[0]); - getStandardChemPotentials(&muss[0]); - getActivities(&actMolal[0]); - - - if (show_thermo) { - sprintf(p, " \n"); - s += p; - sprintf(p, " 1 kg 1 kmol\n"); - s += p; - sprintf(p, " ----------- ------------\n"); - s += p; - sprintf(p, " enthalpy %12.6g %12.4g J\n", - enthalpy_mass(), enthalpy_mole()); - s += p; - sprintf(p, " internal energy %12.6g %12.4g J\n", - intEnergy_mass(), intEnergy_mole()); - s += p; - sprintf(p, " entropy %12.6g %12.4g J/K\n", - entropy_mass(), entropy_mole()); - s += p; - sprintf(p, " Gibbs function %12.6g %12.4g J\n", - gibbs_mass(), gibbs_mole()); - s += p; - sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n", - cp_mass(), cp_mole()); - s += p; - try { - sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n", - cv_mass(), cv_mole()); - s += p; - } - catch(CanteraError) { - sprintf(p, " heat capacity c_v \n"); - s += p; - } - } - - } catch (CanteraError) { - ; - } - return s; - } - } diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h index e9044b1bb..b75594c76 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.h +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -78,18 +78,20 @@ namespace Cantera { *

Specification of Solution Thermodynamic Properties

*
* - * The excess Gibbs free energy + * The excess Gibbs free energy (expressed as an extrinsic thermodynamic + * variable) is given by the following formula: * * \f[ - * G^E = \sum_i \left( H_{Ei} - T S_{Ei} \right) + * G^E = \sum_i \left( H_{Ei} - T S_{Ei} \right) * \f] * \f[ - * H^E_i = X_{Ai} X_{Bi} \left( h_{o,i} + h_{1,i} X_{Bi} \right) + * H^E_i = n X_{Ai} X_{Bi} \left( h_{o,i} + h_{1,i} X_{Bi} \right) * \f] * \f[ - * S^E_i = X_{Ai} X_{Bi} \left( s_{o,i} + s_{1,i} X_{Bi} \right) + * 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. @@ -565,6 +567,18 @@ namespace Cantera { virtual void getPartialMolarEntropies(doublereal* sbar) 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. @@ -579,6 +593,19 @@ namespace Cantera { void getElectrochemPotentials(doublereal* mu) const; + //! Get the array of change in the log activity coefficients with change in state (change temp, change mole fractions) + /*! + * This function is a virtual class, but it first appears in GibbsExcessVPSSTP + * class and derived classes from GibbsExcessVPSSTP. + * + * units = 1/Kelvin + * + * @param dlnActCoeff Output vector of temperature derivatives of the + * log Activity Coefficients. length = m_kk + * + */ + virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeffdT) const; + //! Get the array of temperature derivatives of the log activity coefficients /*! * This function is a virtual class, but it first appears in GibbsExcessVPSSTP @@ -603,16 +630,17 @@ namespace Cantera { * logarithm of the concentration-like variable (i.e. mole fraction, * molality, etc.) 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 + j that concentration-like variable when the derivative of the chemical * potential is taken. * * units = dimensionless * - * @param dlnActCoeffdlnC Output vector of log(mole fraction) + * @param dlnActCoeffdlnX Output vector of log(mole fraction) * derivatives of the log Activity Coefficients. * length = m_kk */ - virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const; + virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const; + virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const; //@} @@ -703,15 +731,6 @@ namespace Cantera { */ void initThermoXML(XML_Node& phaseNode, std::string id); - - - //! returns a summary of the state of the phase as a string - /*! - * @param show_thermo If true, extra information is printed out - * about the thermodynamic state of the system. - */ - virtual std::string report(bool show_thermo = true) const; - private: @@ -761,7 +780,16 @@ namespace Cantera { * derivative of the natural logarithm of the activity coefficients * wrt logarithm of the mole fractions. */ - void s_update_dlnActCoeff_dlnC() const; + void s_update_dlnActCoeff_dlnX() const; + + //! Update the derivative of the log of the activity coefficients + //! wrt log(moles) + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the moles. + */ + void s_update_dlnActCoeff_dlnN() const; private: @@ -802,6 +830,30 @@ namespace Cantera { //! 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 /*! diff --git a/Cantera/src/thermo/MolalityVPSSTP.cpp b/Cantera/src/thermo/MolalityVPSSTP.cpp index 285289a9e..da2a30e60 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.cpp +++ b/Cantera/src/thermo/MolalityVPSSTP.cpp @@ -24,6 +24,7 @@ #include "MolalityVPSSTP.h" +#include using namespace std; namespace Cantera { @@ -910,6 +911,150 @@ namespace Cantera { return s; } + /* + * Format a summary of the mixture state for output. + */ + void MolalityVPSSTP::reportCSV(std::ofstream& csvFile) const { + + + csvFile.precision(3); + int tabS = 15; + int tabM = 30; + int tabL = 40; + try { + if (name() != "") { + csvFile << "\n"+name()+"\n\n"; + } + csvFile << setw(tabL) << "temperature (K) =" << setw(tabS) << temperature() << endl; + csvFile << setw(tabL) << "pressure (Pa) =" << setw(tabS) << pressure() << endl; + csvFile << setw(tabL) << "density (kg/m^3) =" << setw(tabS) << density() << endl; + csvFile << setw(tabL) << "mean mol. weight (amu) =" << setw(tabS) << meanMolecularWeight() << endl; + csvFile << setw(tabL) << "potential (V) =" << setw(tabS) << electricPotential() << endl; + csvFile << endl; + + csvFile << setw(tabL) << "enthalpy (J/kg) = " << setw(tabS) << enthalpy_mass() << setw(tabL) << "enthalpy (J/kmol) = " << setw(tabS) << enthalpy_mole() << endl; + csvFile << setw(tabL) << "internal E (J/kg) = " << setw(tabS) << intEnergy_mass() << setw(tabL) << "internal E (J/kmol) = " << setw(tabS) << intEnergy_mole() << endl; + csvFile << setw(tabL) << "entropy (J/kg) = " << setw(tabS) << entropy_mass() << setw(tabL) << "entropy (J/kmol) = " << setw(tabS) << entropy_mole() << endl; + csvFile << setw(tabL) << "Gibbs (J/kg) = " << setw(tabS) << gibbs_mass() << setw(tabL) << "Gibbs (J/kmol) = " << setw(tabS) << gibbs_mole() << endl; + csvFile << setw(tabL) << "heat capacity c_p (J/K/kg) = " << setw(tabS) << cp_mass() << setw(tabL) << "heat capacity c_p (J/K/kmol) = " << setw(tabS) << cp_mole() << endl; + csvFile << setw(tabL) << "heat capacity c_v (J/K/kg) = " << setw(tabS) << cv_mass() << setw(tabL) << "heat capacity c_v (J/K/kmol) = " << setw(tabS) << cv_mole() << endl; + + csvFile.precision(8); + + int kk = nSpecies(); + double x[kk]; + double molal[kk]; + double mu[kk]; + double muss[kk]; + double aMolal[kk]; + double acMolal[kk]; + double hbar[kk]; + double sbar[kk]; + double ubar[kk]; + double cpbar[kk]; + double vbar[kk]; + vector pNames; + vector data; + + getMoleFractions(x); + pNames.push_back("X"); + data.push_back(x); + try{ + getMolalities(molal); + pNames.push_back("Molal"); + data.push_back(molal); + } + catch (CanteraError) {;} + try{ + getChemPotentials(mu); + pNames.push_back("Chem. Pot. (J/kmol)"); + data.push_back(mu); + } + catch (CanteraError) {;} + try{ + getStandardChemPotentials(muss); + pNames.push_back("Chem. Pot. SS (J/kmol)"); + data.push_back(muss); + } + catch (CanteraError) {;} + try{ + getMolalityActivityCoefficients(acMolal); + pNames.push_back("Molal Act. Coeff."); + data.push_back(acMolal); + } + catch (CanteraError) {;} + try{ + getActivities(aMolal); + pNames.push_back("Molal Activity"); + data.push_back(aMolal); + int iHp = speciesIndex("H+"); + if (iHp >= 0) { + double pH = -log(aMolal[iHp]) / log(10.0); + csvFile << setw(tabL) << "pH = " << setw(tabS) << pH << endl; + } + } + catch (CanteraError) {;} + try{ + getPartialMolarEnthalpies(hbar); + pNames.push_back("Part. Mol Enthalpy (J/kmol)"); + data.push_back(hbar); + } + catch (CanteraError) {;} + try{ + getPartialMolarEntropies(sbar); + pNames.push_back("Part. Mol. Entropy (J/K/kmol)"); + data.push_back(sbar); + } + catch (CanteraError) {;} + try{ + getPartialMolarIntEnergies(ubar); + pNames.push_back("Part. Mol. Energy (J/kmol)"); + data.push_back(ubar); + } + catch (CanteraError) {;} + try{ + getPartialMolarCp(cpbar); + pNames.push_back("Part. Mol. Cp (J/K/kmol"); + data.push_back(cpbar); + } + catch (CanteraError) {;} + try{ + getPartialMolarVolumes(vbar); + pNames.push_back("Part. Mol. Cv (J/K/kmol)"); + data.push_back(vbar); + } + catch (CanteraError) {;} + + csvFile << endl << setw(tabS) << "Species,"; + for ( int i = 0; i < (int)pNames.size(); i++ ){ + csvFile << setw(tabM) << pNames[i] << ","; + } + csvFile << endl; + /* + csvFile.fill('-'); + csvFile << setw(tabS+(tabM+1)*pNames.size()) << "-\n"; + csvFile.fill(' '); + */ + for (int k = 0; k < kk; k++) { + csvFile << setw(tabS) << speciesName(k) + ","; + if (x[k] > SmallNumber) { + for ( int i = 0; i < (int)pNames.size(); i++ ){ + csvFile << setw(tabM) << data[i][k] << ","; + } + csvFile << endl; + } + else{ + for ( int i = 0; i < (int)pNames.size(); i++ ){ + csvFile << setw(tabM) << 0 << ","; + } + csvFile << endl; + } + } + } + catch (CanteraError) { + ; + } + } } diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h index f9eb439d9..60ea18feb 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.h +++ b/Cantera/src/thermo/MolalityVPSSTP.h @@ -800,6 +800,14 @@ namespace Cantera { */ virtual std::string report(bool show_thermo = true) const; + //! returns a summary of the state of the phase to specified + //! comma separated files + /*! + * @param csvFile ofstream file to print comma separated data for + * the phase + */ + virtual void reportCSV(std::ofstream& csvFile) const; + protected: //! Get the array of unscaled non-dimensional molality based diff --git a/Cantera/src/thermo/PDSS_SSVol.cpp b/Cantera/src/thermo/PDSS_SSVol.cpp index 71ae96fd0..4ca9cd46f 100644 --- a/Cantera/src/thermo/PDSS_SSVol.cpp +++ b/Cantera/src/thermo/PDSS_SSVol.cpp @@ -134,14 +134,14 @@ namespace Cantera { m_constMolarVolume = getFloat(*ss, "molarVolume", "toSI"); } else if (model == "temperature_polynomial") { volumeModel_ = cSSVOLUME_TPOLY; - int num = getFloatArray(*ss, TCoeff_, true, "", "volumeTemperaturePolynomial"); + int num = getFloatArray(*ss, TCoeff_, true, "toSI", "volumeTemperaturePolynomial"); if (num != 4) { throw CanteraError("PDSS_SSVol::constructPDSSXML", " Didn't get 4 density polynomial numbers for species " + speciesNode.name()); } } else if (model == "density_temperature_polynomial") { volumeModel_ = cSSVOLUME_DENSITY_TPOLY; - int num = getFloatArray(*ss, TCoeff_, true, "", "densityTemperaturePolynomial"); + int num = getFloatArray(*ss, TCoeff_, true, "toSI", "densityTemperaturePolynomial"); if (num != 4) { throw CanteraError("PDSS_SSVol::constructPDSSXML", " Didn't get 4 density polynomial numbers for species " + speciesNode.name()); diff --git a/Cantera/src/thermo/PureFluidPhase.cpp b/Cantera/src/thermo/PureFluidPhase.cpp index a0748e18c..b5e8b96ff 100644 --- a/Cantera/src/thermo/PureFluidPhase.cpp +++ b/Cantera/src/thermo/PureFluidPhase.cpp @@ -17,6 +17,7 @@ #include "../../../ext/tpx/utils.h" #include +#include namespace Cantera { @@ -427,8 +428,144 @@ namespace Cantera { } return s; } + + /* + * Format a summary of the mixture state for output. + */ + void PureFluidPhase::reportCSV(std::ofstream& csvFile) const { + + + csvFile.precision(3); + int tabS = 15; + int tabM = 30; + int tabL = 40; + try { + if (name() != "") { + csvFile << "\n"+name()+"\n\n"; + } + csvFile << setw(tabL) << "temperature (K) =" << setw(tabS) << temperature() << endl; + csvFile << setw(tabL) << "pressure (Pa) =" << setw(tabS) << pressure() << endl; + csvFile << setw(tabL) << "density (kg/m^3) =" << setw(tabS) << density() << endl; + csvFile << setw(tabL) << "mean mol. weight (amu) =" << setw(tabS) << meanMolecularWeight() << endl; + csvFile << setw(tabL) << "potential (V) =" << setw(tabS) << electricPotential() << endl; + if (eosType() == cPureFluid) { + double xx = ((PureFluidPhase *) (this))->vaporFraction(); + csvFile << setw(tabL) << "vapor fraction = " << setw(tabS) << xx << endl; + } + csvFile << endl; + + csvFile << setw(tabL) << "enthalpy (J/kg) = " << setw(tabS) << enthalpy_mass() << setw(tabL) << "enthalpy (J/kmol) = " << setw(tabS) << enthalpy_mole() << endl; + csvFile << setw(tabL) << "internal E (J/kg) = " << setw(tabS) << intEnergy_mass() << setw(tabL) << "internal E (J/kmol) = " << setw(tabS) << intEnergy_mole() << endl; + csvFile << setw(tabL) << "entropy (J/kg) = " << setw(tabS) << entropy_mass() << setw(tabL) << "entropy (J/kmol) = " << setw(tabS) << entropy_mole() << endl; + csvFile << setw(tabL) << "Gibbs (J/kg) = " << setw(tabS) << gibbs_mass() << setw(tabL) << "Gibbs (J/kmol) = " << setw(tabS) << gibbs_mole() << endl; + csvFile << setw(tabL) << "heat capacity c_p (J/K/kg) = " << setw(tabS) << cp_mass() << setw(tabL) << "heat capacity c_p (J/K/kmol) = " << setw(tabS) << cp_mole() << endl; + csvFile << setw(tabL) << "heat capacity c_v (J/K/kg) = " << setw(tabS) << cv_mass() << setw(tabL) << "heat capacity c_v (J/K/kmol) = " << setw(tabS) << cv_mole() << endl; + + csvFile.precision(8); + + int kk = nSpecies(); + std::vector x(kk, 0.0); + std::vector y(kk, 0.0); + std::vector mu(kk, 0.0); + std::vector a(kk, 0.0); + std::vector ac(kk, 0.0); + std::vector hbar(kk, 0.0); + std::vector sbar(kk, 0.0); + std::vector ubar(kk, 0.0); + std::vector cpbar(kk, 0.0); + std::vector vbar(kk, 0.0); + vector pNames; + vector data; + + getMoleFractions(DATA_PTR(x)); + pNames.push_back("X"); + data.push_back(DATA_PTR(x)); + try{ + getMassFractions(DATA_PTR(y)); + pNames.push_back("Y"); + data.push_back(DATA_PTR(y)); + } + catch (CanteraError) {;} + try{ + getChemPotentials(DATA_PTR(mu)); + pNames.push_back("Chem. Pot (J/kmol)"); + data.push_back(DATA_PTR(mu)); + } + catch (CanteraError) {;} + try{ + getActivities(DATA_PTR(a)); + pNames.push_back("Activity"); + data.push_back(DATA_PTR(a)); + } + catch (CanteraError) {;} + try{ + getActivityCoefficients(DATA_PTR(ac)); + pNames.push_back("Act. Coeff."); + data.push_back(DATA_PTR(ac)); + } + catch (CanteraError) {;} + try{ + getPartialMolarEnthalpies(DATA_PTR(hbar)); + pNames.push_back("Part. Mol Enthalpy (J/kmol)"); + data.push_back(DATA_PTR(hbar)); + } + catch (CanteraError) {;} + try{ + getPartialMolarEntropies(DATA_PTR(sbar)); + pNames.push_back("Part. Mol. Entropy (J/K/kmol)"); + data.push_back(DATA_PTR(sbar)); + } + catch (CanteraError) {;} + try{ + getPartialMolarIntEnergies(DATA_PTR(ubar)); + pNames.push_back("Part. Mol. Energy (J/kmol)"); + data.push_back(DATA_PTR(ubar)); + } + catch (CanteraError) {;} + try{ + getPartialMolarCp(DATA_PTR(cpbar)); + pNames.push_back("Part. Mol. Cp (J/K/kmol"); + data.push_back(DATA_PTR(cpbar)); + } + catch (CanteraError) {;} + try{ + getPartialMolarVolumes(DATA_PTR(vbar)); + pNames.push_back("Part. Mol. Cv (J/K/kmol)"); + data.push_back(DATA_PTR(vbar)); + } + catch (CanteraError) {;} + + csvFile << endl << setw(tabS) << "Species,"; + for ( int i = 0; i < (int)pNames.size(); i++ ){ + csvFile << setw(tabM) << pNames[i] << ","; + } + csvFile << endl; + /* + csvFile.fill('-'); + csvFile << setw(tabS+(tabM+1)*pNames.size()) << "-\n"; + csvFile.fill(' '); + */ + for (int k = 0; k < kk; k++) { + csvFile << setw(tabS) << speciesName(k) + ","; + if (x[k] > SmallNumber) { + for ( int i = 0; i < (int)pNames.size(); i++ ){ + csvFile << setw(tabM) << data[i][k] << ","; + } + csvFile << endl; + } + else{ + for ( int i = 0; i < (int)pNames.size(); i++ ){ + csvFile << setw(tabM) << 0 << ","; + } + csvFile << endl; + } + } + } + catch (CanteraError) { + ; + } + } - } #endif // WITH_PURE_FLUIDS diff --git a/Cantera/src/thermo/PureFluidPhase.h b/Cantera/src/thermo/PureFluidPhase.h index 2e2c626c1..74626fe22 100644 --- a/Cantera/src/thermo/PureFluidPhase.h +++ b/Cantera/src/thermo/PureFluidPhase.h @@ -303,6 +303,14 @@ namespace Cantera { */ virtual std::string report(bool show_thermo = true) const; + //! returns a summary of the state of the phase to specified + //! comma separated files + /*! + * @param csvFile ofstream file to print comma separated data for + * the phase + */ + virtual void reportCSV(std::ofstream& csvFile) const; + protected: //! Main call to the tpx level to set the state of the system diff --git a/Cantera/src/thermo/State.cpp b/Cantera/src/thermo/State.cpp index 6eeec8ce8..a25d53701 100644 --- a/Cantera/src/thermo/State.cpp +++ b/Cantera/src/thermo/State.cpp @@ -193,6 +193,10 @@ namespace Cantera { return density()/meanMolecularWeight(); } + doublereal State::molarVolume() const { + return 1.0/molarDensity(); + } + void State::setConcentrations(const doublereal* const conc) { int k; doublereal sum = 0.0, norm = 0.0; diff --git a/Cantera/src/thermo/State.h b/Cantera/src/thermo/State.h index d05941add..ca2826df9 100644 --- a/Cantera/src/thermo/State.h +++ b/Cantera/src/thermo/State.h @@ -318,6 +318,9 @@ namespace Cantera { /// Molar density (kmol/m^3). doublereal molarDensity() const; + /// Molar density (kmol/m^3). + doublereal molarVolume() const; + //! Set the internally storred density (kg/m^3) of the phase /*! * Note the density of a phase is an indepedent variable. diff --git a/Cantera/src/thermo/ThermoPhase.cpp b/Cantera/src/thermo/ThermoPhase.cpp index eb0e60250..afc819617 100644 --- a/Cantera/src/thermo/ThermoPhase.cpp +++ b/Cantera/src/thermo/ThermoPhase.cpp @@ -21,6 +21,7 @@ #endif #include "ThermoPhase.h" +#include //@{ #ifndef MAX @@ -49,11 +50,13 @@ namespace Cantera { ThermoPhase::~ThermoPhase() { for (int k = 0; k < m_kk; k++) { - if (!m_speciesData[k]) { + if (m_speciesData[k]) { delete m_speciesData[k]; + m_speciesData[k] = 0; } } delete m_spthermo; + m_spthermo = 0; } /** @@ -95,8 +98,9 @@ namespace Cantera { * We need to destruct first */ for (int k = 0; k < m_kk; k++) { - if (!m_speciesData[k]) { + if (m_speciesData[k]) { delete m_speciesData[k]; + m_speciesData[k] = 0; } } if (m_spthermo) { @@ -1132,5 +1136,153 @@ namespace Cantera { return s; } + /* + * Format a summary of the mixture state for output. + */ + void ThermoPhase::reportCSV(std::ofstream& csvFile) const { + + csvFile.precision(3); + int tabS = 15; + int tabM = 30; + int tabL = 40; + try { + if (name() != "") { + csvFile << "\n"+name()+"\n\n"; + } + csvFile << setw(tabL) << "temperature (K) =" << setw(tabS) << temperature() << endl; + csvFile << setw(tabL) << "pressure (Pa) =" << setw(tabS) << pressure() << endl; + csvFile << setw(tabL) << "density (kg/m^3) =" << setw(tabS) << density() << endl; + csvFile << setw(tabL) << "mean mol. weight (amu) =" << setw(tabS) << meanMolecularWeight() << endl; + csvFile << setw(tabL) << "potential (V) =" << setw(tabS) << electricPotential() << endl; + csvFile << endl; + + csvFile << setw(tabL) << "enthalpy (J/kg) = " << setw(tabS) << enthalpy_mass() << setw(tabL) + << "enthalpy (J/kmol) = " << setw(tabS) << enthalpy_mole() << endl; + csvFile << setw(tabL) << "internal E (J/kg) = " << setw(tabS) << intEnergy_mass() << setw(tabL) + << "internal E (J/kmol) = " << setw(tabS) << intEnergy_mole() << endl; + csvFile << setw(tabL) << "entropy (J/kg) = " << setw(tabS) << entropy_mass() << setw(tabL) + << "entropy (J/kmol) = " << setw(tabS) << entropy_mole() << endl; + csvFile << setw(tabL) << "Gibbs (J/kg) = " << setw(tabS) << gibbs_mass() << setw(tabL) + << "Gibbs (J/kmol) = " << setw(tabS) << gibbs_mole() << endl; + csvFile << setw(tabL) << "heat capacity c_p (J/K/kg) = " << setw(tabS) << cp_mass() + << setw(tabL) << "heat capacity c_p (J/K/kmol) = " << setw(tabS) << cp_mole() << endl; + csvFile << setw(tabL) << "heat capacity c_v (J/K/kg) = " << setw(tabS) << cv_mass() + << setw(tabL) << "heat capacity c_v (J/K/kmol) = " << setw(tabS) << cv_mole() << endl; + + csvFile.precision(8); + + int kk = nSpecies(); + doublereal *x = new doublereal[kk]; + doublereal *y = new doublereal[kk]; + doublereal *mu = new doublereal[kk]; + doublereal *a = new doublereal[kk]; + doublereal *ac = new doublereal[kk]; + doublereal *hbar = new doublereal[kk]; + doublereal *sbar = new doublereal[kk]; + doublereal *ubar = new doublereal[kk]; + doublereal *cpbar= new doublereal[kk]; + doublereal *vbar = new doublereal[kk]; + std::vector pNames; + std::vector data; + + getMoleFractions(x); + pNames.push_back("X"); + data.push_back(x); + try{ + getMassFractions(y); + pNames.push_back("Y"); + data.push_back(y); + } + catch (CanteraError) {;} + try{ + getChemPotentials(mu); + pNames.push_back("Chem. Pot (J/kmol)"); + data.push_back(mu); + } + catch (CanteraError) {;} + try{ + getActivities(a); + pNames.push_back("Activity"); + data.push_back(a); + } + catch (CanteraError) {;} + try{ + getActivityCoefficients(ac); + pNames.push_back("Act. Coeff."); + data.push_back(ac); + } + catch (CanteraError) {;} + try{ + getPartialMolarEnthalpies(hbar); + pNames.push_back("Part. Mol Enthalpy (J/kmol)"); + data.push_back(hbar); + } + catch (CanteraError) {;} + try{ + getPartialMolarEntropies(sbar); + pNames.push_back("Part. Mol. Entropy (J/K/kmol)"); + data.push_back(sbar); + } + catch (CanteraError) {;} + try{ + getPartialMolarIntEnergies(ubar); + pNames.push_back("Part. Mol. Energy (J/kmol)"); + data.push_back(ubar); + } + catch (CanteraError) {;} + try{ + getPartialMolarCp(cpbar); + pNames.push_back("Part. Mol. Cp (J/K/kmol"); + data.push_back(cpbar); + } + catch (CanteraError) {;} + try{ + getPartialMolarVolumes(vbar); + pNames.push_back("Part. Mol. Cv (J/K/kmol)"); + data.push_back(vbar); + } + catch (CanteraError) {;} + + csvFile << endl << setw(tabS) << "Species,"; + for ( int i = 0; i < (int)pNames.size(); i++ ){ + csvFile << setw(tabM) << pNames[i] << ","; + } + csvFile << endl; + /* + csvFile.fill('-'); + csvFile << setw(tabS+(tabM+1)*pNames.size()) << "-\n"; + csvFile.fill(' '); + */ + for (int k = 0; k < kk; k++) { + csvFile << setw(tabS) << speciesName(k) + ","; + if (x[k] > SmallNumber) { + for (int i = 0; i < (int)pNames.size(); i++) { + csvFile << setw(tabM) << data[i][k] << ","; + } + csvFile << endl; + } else { + for (int i = 0; i < (int)pNames.size(); i++) { + csvFile << setw(tabM) << 0 << ","; + } + csvFile << endl; + } + } + delete [] x; + delete [] y; + delete [] mu; + delete [] a; + delete [] ac; + delete [] hbar; + delete [] sbar; + delete [] ubar; + delete [] cpbar; + delete [] vbar; + + } + catch (CanteraError) { + ; + } + } + } diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h index 3d5ddb577..c52e74e75 100644 --- a/Cantera/src/thermo/ThermoPhase.h +++ b/Cantera/src/thermo/ThermoPhase.h @@ -883,6 +883,21 @@ namespace Cantera { } + //! Get the change in activity coefficients w.r.t. change in state + //! (temp, mole fraction, etc.) + /*! + * This function is a virtual method. For ideal mixtures + * (unity activity coefficients), this can gradX/X. + * + * @param dT Input of temperature change + * @param dX Input vector of changes in mole fraction. length = m_kk + * @param dlnActCoeff Output vector of derivatives of the + * log Activity Coefficients. length = m_kk + */ + virtual void getdlnActCoeff(const doublereal dT, const doublereal * const dX, doublereal *dlnActCoeff) const { + err("getdlnActCoeff"); + } + //! Get the array of log concentration-like derivatives of the //! log activity coefficients /*! @@ -890,19 +905,41 @@ namespace Cantera { * (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. + * logarithm of the concentration-like variable (i.e. mole fraction) + * 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 dlnActCoeffdlnC Output vector of derivatives of the + * @param dlnActCoeffdlnX Output vector of derivatives of the * log Activity Coefficients. length = m_kk */ - virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const { - err("getdlnActCoeffdlnC"); + virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const { + err("getdlnActCoeffdlnX"); + } + + //! 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 Output vector of derivatives of the + * log Activity Coefficients. length = m_kk + */ + virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const { + err("getdlnActCoeffdlnN"); } @@ -2074,6 +2111,13 @@ namespace Cantera { * about the thermodynamic state of the system. */ virtual std::string report(bool show_thermo = true) const; + + //! returns a summary of the state of the phase to a comma separated file + /*! + * @param csvFile ofstream file to print comma separated data for + * the phase + */ + virtual void reportCSV(std::ofstream& csvFile) const; protected: diff --git a/Cantera/src/thermo/VPSSMgr.cpp b/Cantera/src/thermo/VPSSMgr.cpp index 64ec7e22c..1204bce2c 100644 --- a/Cantera/src/thermo/VPSSMgr.cpp +++ b/Cantera/src/thermo/VPSSMgr.cpp @@ -280,7 +280,8 @@ namespace Cantera { void VPSSMgr::getStandardVolumes_ref(doublereal *vol) const{ - err("getStandardVolumes_ref"); + getStandardVolumes(vol); + //err("getStandardVolumes_ref"); } /*****************************************************************/ @@ -359,6 +360,7 @@ namespace Cantera { m_sss_R.resize(m_kk, 0.0); m_Vss.resize(m_kk, 0.0); + // Storage used by the PDSS objects to store their // answers. mPDSS_h0_RT.resize(m_kk, 0.0); diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h index c13412fe4..7990b2cff 100644 --- a/Cantera/src/thermo/VPStandardStateTP.h +++ b/Cantera/src/thermo/VPStandardStateTP.h @@ -127,19 +127,41 @@ namespace Cantera { * (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. + * 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 dlnActCoeffdlnC Output vector of derivatives of the + * @param dlnActCoeffdlnN Output vector of derivatives of the * log Activity Coefficients. length = m_kk */ - virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const { - err("getdlnActCoeffdlnC"); + virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const { + err("getdlnActCoeffdlnN"); + } + + //! 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. mole fraction) + * 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 dlnActCoeffdlnX Output vector of derivatives of the + * log Activity Coefficients. length = m_kk + */ + virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const { + err("getdlnActCoeffdlnX"); } diff --git a/Cantera/src/thermo/WaterProps.cpp b/Cantera/src/thermo/WaterProps.cpp index 9ba777163..f9a4421a0 100644 --- a/Cantera/src/thermo/WaterProps.cpp +++ b/Cantera/src/thermo/WaterProps.cpp @@ -501,6 +501,7 @@ namespace Cantera { 0.978197, 0.579829, -0.202354}; + //! parameter const doublereal Hij[6][7] = {