research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Worked on MolarityIonicVPSSTP, filling it in a little bit.

Browse source 
Changed GibbsExcessVPSSTP so that getLnActivityCoefficient() is central.
The straight activity coefficients may not be representable within
machine limits.
This commit is contained in:
Harry Moffat 2011-11-01 21:09:58 +00:00
parent 5d1e55596d
commit 5f01cc30bb
8 changed files with 713 additions and 38 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
Cantera/src/thermo/GibbsExcessVPSSTP.cpp
View file

@ -233,7 +233,29 @@ namespace Cantera {
}
}
//====================================================================================================================
// 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 GibbsExcessVPSSTP::getActivityCoefficients(doublereal * const ac) const {
getLnActivityCoefficients(ac);
// Protect against roundoff when taking exponentials
for (int k = 0; k < m_kk; k++) {
if (ac[k] > 700.) {
ac[k] = exp(700.);
} else if (ac[k] < -700.) {
ac[k] = exp(-700);
} else {
ac[k] = exp(ac[k]);
}
}
}
//====================================================================================================================
void GibbsExcessVPSSTP::getElectrochemPotentials(doublereal* mu) const {
getChemPotentials(mu);
double ve = Faraday * electricPotential();

7
Cantera/src/thermo/GibbsExcessVPSSTP.h
View file

@ -289,6 +289,13 @@ namespace Cantera {
*/
virtual void getActivities(doublereal* ac) const;
//! Get the array of non-dimensional molar-based ln activity coefficients at
//! the current solution temperature, pressure, and solution concentration.
/*!
* @param lnac Output vector of ln activity coefficients. Length: m_kk.
*/
virtual void getActivityCoefficients(doublereal * const ac) const;
//! Get the array of temperature derivatives of the log activity coefficients
/*!

10
Cantera/src/thermo/MargulesVPSSTP.cpp
View file

@ -376,12 +376,12 @@ namespace Cantera {
return 0.0;
}
//====================================================================================================================
// Get the array of non-dimensional molar-based activity coefficients at
// Get the array of non-dimensional molar-based ln activity coefficients at
// the current solution temperature, pressure, and solution concentration.
/*
* @param ac Output vector of activity coefficients. Length: m_kk.
* @param lnac Output vector of ln activity coefficients. Length: m_kk.
*/
void MargulesVPSSTP::getActivityCoefficients(doublereal* ac) const {
void MargulesVPSSTP::getLnActivityCoefficients(doublereal* lnac) const {
/*
* Update the activity coefficients
*/
@ -391,10 +391,10 @@ namespace Cantera {
* take the exp of the internally storred coefficients.
*/
for (int k = 0; k < m_kk; k++) {
ac[k] = exp(lnActCoeff_Scaled_[k]);
lnac[k] = lnActCoeff_Scaled_[k];
}
}
//====================================================================================================================
/*
* ------------ Partial Molar Properties of the Solution ------------
*/

8
Cantera/src/thermo/MargulesVPSSTP.h
View file

@ -536,15 +536,13 @@ namespace Cantera {
*/
virtual doublereal logStandardConc(int k=0) const;
//! Get the array of non-dimensional molar-based activity coefficients at
//! Get the array of non-dimensional molar-based ln activity coefficients at
//! the current solution temperature, pressure, and solution concentration.
/*!
* @param ac Output vector of activity coefficients. Length: m_kk.
* @param ac Output vector of ln activity coefficients. Length: m_kk.
*/
virtual void getActivityCoefficients(doublereal* ac) const;
virtual void getLnActivityCoefficients(doublereal* lnac) const;
//@}
/// @name Partial Molar Properties of the Solution

453
Cantera/src/thermo/MolarityIonicVPSSTP.cpp
View file

@ -23,12 +23,13 @@
#include "MolarityIonicVPSSTP.h"
#include "ThermoFactory.h"
#include <cmath>
using namespace std;
namespace Cantera {
static const double xxSmall = 1.0E-150;
//====================================================================================================================
/*
* Default constructor.
@ -44,6 +45,40 @@ namespace Cantera {
numPassThroughSpecies_(0),
neutralPBindexStart(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.
*/
MolarityIonicVPSSTP::MolarityIonicVPSSTP(std::string inputFile, std::string id) :
GibbsExcessVPSSTP(),
PBType_(PBTYPE_PASSTHROUGH),
numPBSpecies_(m_kk),
indexSpecialSpecies_(-1),
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralPBindexStart(0)
{
constructPhaseFile(inputFile, id);
}
//====================================================================================================================
MolarityIonicVPSSTP::MolarityIonicVPSSTP(XML_Node& phaseRoot, std::string id) :
GibbsExcessVPSSTP(),
PBType_(PBTYPE_PASSTHROUGH),
numPBSpecies_(m_kk),
indexSpecialSpecies_(-1),
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralPBindexStart(0)
{
constructPhaseXML(phaseRoot, id);
}
//====================================================================================================================
/*
@ -130,8 +165,120 @@ namespace Cantera {
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 MolarityIonicVPSSTP::constructPhaseFile(std::string inputFile, std::string id) {
if ((int) inputFile.size() == 0) {
throw CanteraError("MolarityIonicVPSSTP:constructPhaseFile",
"input file is null");
}
string path = findInputFile(inputFile);
std::ifstream fin(path.c_str());
if (!fin) {
throw CanteraError("MolarityIonicVPSSTP: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("MolarityIonicVPSSTP: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 MolarityIonicVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id) {
string stemp;
if ((int) id.size() > 0) {
string idp = phaseNode.id();
if (idp != id) {
throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML",
"phasenode and Id are incompatible");
}
}
/*
* Find the Thermo XML node
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
/*
* Make sure that the thermo model is MolarityIonic
*/
stemp = thermoNode.attrib("model");
string formString = lowercase(stemp);
if (formString != "molarityionicvpss") {
throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML",
"model name isn't MolarityIonicVPSS: " + 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("MolarityIonicVPSSTP::constructPhaseXML","importPhase failed ");
}
}
//====================================================================================================================
/*
* ------------ Molar Thermodynamic Properties ----------------------
*/
@ -140,16 +287,78 @@ namespace Cantera {
/*
* - 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 MolarityIonicVPSSTP::getActivityConcentrations(doublereal* c) const {
getActivities(c);
}
//====================================================================================================================
doublereal MolarityIonicVPSSTP::standardConcentration(int k) const {
err("standardConcentration");
return -1.0;
}
//====================================================================================================================
doublereal MolarityIonicVPSSTP::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 MolarityIonicVPSSTP::getLnActivityCoefficients(doublereal* lnac) const {
/*
* Update the activity coefficients
*/
s_update_lnActCoeff();
/*
* take the exp of the internally storred coefficients.
*/
for (int k = 0; k < m_kk; k++) {
lnac[k] = lnActCoeff_Scaled_[k];
}
}
//====================================================================================================================
void MolarityIonicVPSSTP::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]);
}
}
//====================================================================================================================
@ -160,6 +369,145 @@ namespace Cantera {
mu[k] += ve*charge(k);
}
}
//====================================================================================================================
// 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 MolarityIonicVPSSTP::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 MolarityIonicVPSSTP::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 MolarityIonicVPSSTP::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;
}
}
// 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 MolarityIonicVPSSTP::getPartialMolarVolumes(doublereal* vbar) const {
int iK;
/*
* Get the standard state values in m^3 kmol-1
*/
getStandardVolumes(vbar);
for ( iK = 0; iK < m_kk; iK++ ){
vbar[iK] += 0.0;
}
}
//====================================================================================================================
void MolarityIonicVPSSTP::calcPseudoBinaryMoleFractions() const {
int k;
@ -241,6 +589,38 @@ namespace Cantera {
}
}
//====================================================================================================================
// Update the activity coefficients
/*
* This function will be called to update the internally storred
* natural logarithm of the activity coefficients
*
*/
void MolarityIonicVPSSTP::s_update_lnActCoeff() const {
int k;
for (k = 0; k < m_kk; k++) {
lnActCoeff_Scaled_[k] = 0.0;
}
}
//====================================================================================================================
void MolarityIonicVPSSTP::s_update_dlnActCoeff_dT() const {
}
//====================================================================================================================
// Internal routine that calculates the derivative of the activity coefficients wrt
// the mole fractions.
/*
* This routine calculates the the derivative of the activity coefficients wrt to mole fraction
* with all other mole fractions held constant. This is strictly not permitted. However, if the
* resulting matrix is multiplied by a permissible deltaX vector then everything is ok.
*
* This is the natural way to handle concentration derivatives in this routine.
*/
void MolarityIonicVPSSTP::s_update_dlnActCoeff_dX_() const {
}
//====================================================================================================================
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
@ -323,11 +703,70 @@ namespace Cantera {
* with the correct id.
*/
void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id) {
std::string subname = "MolarityIonicVPSSTP::initThermoXML";
std::string stemp;
/*
* Check on the thermo field. Must have:
* <thermo model="MolarityIonic" />
*/
XML_Node& thermoNode = phaseNode.child("thermo");
std::string mStringa = thermoNode.attrib("model");
std::string mString = lowercase(mStringa);
if (mString != "molarityionicvpss") {
throw CanteraError(subname.c_str(),
"Unknown thermo model: " + mStringa + " - This object only knows \"MolarityIonicVPSS\" ");
}
/*
* 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;
std::string mStringa = acNode.attrib("model");
std::string mString = lowercase(mStringa);
// if (mString != "redlich-kister") {
// 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();
std::string nodeName = lowercase(stemp);
/*
* Process a binary interaction
*/
if (nodeName == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
}
/*
* Go down the chain
*/
GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
}
//====================================================================================================================
// Process an XML node called "binaryNeutralSpeciesParameters"
/*
* This node contains all of the parameters necessary to describe
* a single binary interaction. This function reads the XML file and writes the coefficients
* it finds to an internal data structures.
*/
void MolarityIonicVPSSTP::readXMLBinarySpecies(XML_Node &xmLBinarySpecies) {
std::string xname = xmLBinarySpecies.name();
}
//====================================================================================================================
/*
* Format a summary of the mixture state for output.
*/

222
Cantera/src/thermo/MolarityIonicVPSSTP.h
View file

@ -38,8 +38,7 @@ namespace Cantera {
* thermodynamic properties that are further based on
* expressing the Excess Gibbs free energy as a function of
* the mole fractions (or pseudo mole fractions) of the consitituents.
* This category is the workhorse for describing ionic systems which
* are not on the molality scale.
* This category is the workhorse for describing ionic systems which are not on the molality scale.
*
* This class adds additional functions onto the %ThermoPhase interface
* that handles the calculation of the excess Gibbs free energy. The %ThermoPhase
@ -56,7 +55,7 @@ namespace Cantera {
* it is expected that there exists a charge balance at all times.
* One of the ions must be a "special ion" in the sense that its' thermodynamic
* functions are set to zero, and the thermo functions of all other
* ions are based on a valuation relative to the special ion.
* ions are based on a valuation relative to that special ion.
*
*/
class MolarityIonicVPSSTP : public GibbsExcessVPSSTP {
@ -74,6 +73,31 @@ namespace Cantera {
*/
MolarityIonicVPSSTP();
//! Construct and initialize a MolarityIonicVPSSTP 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.
*/
MolarityIonicVPSSTP(std::string inputFile, std::string id = "");
//! Construct and initialize a MolarityIonicVPSSTP 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)
*/
MolarityIonicVPSSTP(XML_Node& phaseRef, std::string id = "");
//! Copy constructor
/*!
* Note this stuff will not work until the underlying phase
@ -118,7 +142,46 @@ namespace Cantera {
*/
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);
/**
* @}
@ -165,7 +228,24 @@ namespace Cantera {
* @{
*/
//! 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;
/**
@ -191,6 +271,15 @@ namespace Cantera {
*/
virtual doublereal logStandardConc(int k=0) const;
//! Get the array of non-dimensional molar-based ln activity coefficients at
//! the current solution temperature, pressure, and solution concentration.
/*!
* @param lnac Output vector of ln activity coefficients. Length: m_kk.
*/
virtual void getLnActivityCoefficients(doublereal* ac) const;
@ -198,6 +287,16 @@ namespace Cantera {
/// @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;
/**
* Get the species electrochemical potentials.
@ -212,7 +311,79 @@ namespace Cantera {
*/
void getElectrochemPotentials(doublereal* mu) 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;
//@}
/// @name Properties of the Standard State of the Species in the Solution
//@{
@ -318,11 +489,48 @@ namespace Cantera {
private:
//! Initialize lengths of local variables after all species have
//! been identified.
//! Initialize lengths of local variables after all species have been identified.
void initLengths();
//! Process an XML node called "binaryNeutralSpeciesParameters"
/*!
* This node contains all of the parameters necessary to describe
* the Redlich-Kister 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);
//! 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;
//! Internal routine that calculates the derivative of the activity coefficients wrt
//! the mole fractions.
/*!
* This routine calculates the the derivative of the activity coefficients wrt to mole fraction
* with all other mole fractions held constant. This is strictly not permitted. However, if the
* resulting matrix is multiplied by a permissible deltaX vector then everything is ok.
*
* This is the natural way to handle concentration derivatives in this routine.
*/
void s_update_dlnActCoeff_dX_() const;
private:
//! Error function

20
Cantera/src/thermo/RedlichKisterVPSSTP.cpp
View file

@ -387,7 +387,7 @@ namespace Cantera {
/*
* @param ac Output vector of activity coefficients. Length: m_kk.
*/
void RedlichKisterVPSSTP::getActivityCoefficients(doublereal* ac) const {
void RedlichKisterVPSSTP::getLnActivityCoefficients(doublereal* lnac) const {
/*
* Update the activity coefficients
*/
@ -397,7 +397,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]);
lnac[k] = lnActCoeff_Scaled_[k];
}
}
//====================================================================================================================
@ -728,6 +728,12 @@ namespace Cantera {
doublereal RT = GasConstant * T;
fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk);
/*
* Scaling: I moved the division of RT higher so that we are always dealing with G/RT dimensionless terms
* within the routine. There is a severe problem with roundoff error in these calculations. The
* dimensionless terms help.
*/
for (int i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
@ -744,7 +750,7 @@ namespace Cantera {
doublereal sumMm1 = 0.0;
doublereal sum2 = 0.0;
for (m = 0; m < N; m++) {
doublereal A_ge = he_vec[m] - T * se_vec[m];
doublereal A_ge = (he_vec[m] - T * se_vec[m]) / RT;
sum += A_ge * poly;
sum2 += A_ge * (m + 1) * poly;
poly *= deltaX;
@ -771,7 +777,7 @@ namespace Cantera {
double polyk = 1.0;
double fac = 2.0 * XA - 1.0;
for (m = 0; m < N; m++) {
doublereal A_ge = he_vec[m] - T * se_vec[m];
doublereal A_ge = (he_vec[m] - T * se_vec[m]) / RT;
lnA += A_ge * oneMXA * oneMXA * polyk * (1.0 + 2.0 * XA * m / fac);
lnB += A_ge * XA * XA * polyk * (1.0 - 2.0 * oneMXA * m / fac);
polyk *= fac;
@ -783,9 +789,7 @@ namespace Cantera {
#endif
}
for (k = 0; k < m_kk; k++) {
lnActCoeff_Scaled_[k] /= RT;
}
}
//===================================================================================================================
// Update the derivative of the log of the activity coefficients wrt T
@ -1117,7 +1121,7 @@ namespace Cantera {
double poly1mk = fac;
for (m = 0; m < N; m++) {
doublereal A_ge = he_vec[m];
doublereal A_ge = he_vec[m] - T * se_vec[m];
Volts += A_ge * ( polykp1 - (2.0 * XA * m * (1.0-XA) ) / poly1mk );
polykp1 *= fac;
poly1mk /= fac;

9
Cantera/src/thermo/RedlichKisterVPSSTP.h
View file

@ -533,15 +533,12 @@ namespace Cantera {
*/
virtual doublereal logStandardConc(int k=0) const;
//! Get the array of non-dimensional molar-based activity coefficients at
//! Get the array of non-dimensional molar-based ln activity coefficients at
//! the current solution temperature, pressure, and solution concentration.
/*!
* @param ac Output vector of activity coefficients. Length: m_kk.
* @param lnac Output vector of ln activity coefficients. Length: m_kk.
*/
virtual void getActivityCoefficients(doublereal* ac) const;
virtual void getLnActivityCoefficients(doublereal* ac) const;
//@}
/// @name Partial Molar Properties of the Solution