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

Merged thermo directory with the LiquidTransportDevelop

Browse source 
thermo directory. 
  The biggest change is the addition of the derivative routines
of the activity wrt mole fraction and mole number.
Also There is an addition of reportCSV() routine, which is 
a method to create a comma separated file from ThermoPhase routines.
The rest are small issues.
This commit is contained in:
Harry Moffat 2010-05-09 03:18:33 +00:00
parent 6089ab2d68
commit ef7728b1d3
21 changed files with 1406 additions and 394 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -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;

1
Cantera/src/thermo/Elements.cpp
View file

@ -328,7 +328,6 @@ namespace Cantera {
}
doublereal Elements::entropyElement298(int m) const {
AssertThrowMsg(m_entropy298[m] != ENTROPY298_UNKNOWN,
"Elements::entropy298",

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

@ -23,6 +23,7 @@
#include "GibbsExcessVPSSTP.h"
#include <iomanip>
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 <not implemented> \n");
s += p;
}
}
} catch (CanteraError) {
;
}
return s;
}
}

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

@ -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.
*
* <H3> SetState Strategy </H3>
*
* 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<doublereal> 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<doublereal> dlnActCoeffdlnC_Scaled_;
mutable std::vector<doublereal> 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<doublereal> dlnActCoeffdlnX_Scaled_;
//! Temporary storage space that is fair game
mutable std::vector<doublereal> m_pp;

1
Cantera/src/thermo/IdealGasPhase.h
View file

@ -330,6 +330,7 @@ namespace Cantera {
*/
IdealGasPhase(const IdealGasPhase &right);
//! Asignment operator
/*!
* Assignment operator for the object. Constructed

420
Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp
View file

@ -28,6 +28,7 @@
#include "mix_defs.h"
#include <cmath>
#include <iomanip>
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<int>& 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<std::string>& elnamesVN ,
@ -1136,9 +1299,12 @@ namespace Cantera {
std::vector<double> elemVectorI(nElementsI);
vector<doublereal> 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<GibbsExcessVPSSTP *>(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<GibbsExcessVPSSTP *>(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 <not implemented> \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<GibbsExcessVPSSTP *>(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;
}
}

108
Cantera/src/thermo/IonsFromNeutralVPSSTP.h
View file

@ -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<int>& 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<int>& cation){cation=cationList_;}
virtual void getAnionList(std::vector<int>& anion){anion=anionList_;}
virtual void getSpeciesNames(std::vector<std::string>& 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<doublereal> muNeutralMolecule_;
mutable std::vector<doublereal> gammaNeutralMolecule_;
mutable std::vector<doublereal> dlnActCoeff_NeutralMolecule_;
mutable std::vector<doublereal> dlnActCoeffdT_NeutralMolecule_;
mutable std::vector<doublereal> dlnActCoeffdlnC_NeutralMolecule_;
mutable std::vector<doublereal> dlnActCoeffdlnX_NeutralMolecule_;
mutable std::vector<doublereal> dlnActCoeffdlnN_NeutralMolecule_;
};

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

@ -19,7 +19,7 @@
#include "MargulesVPSSTP.h"
#include "ThermoFactory.h"
#include <iomanip>
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 <not implemented> \n");
s += p;
}
}
} catch (CanteraError) {
;
}
return s;
}
}

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

@ -78,18 +78,20 @@ namespace Cantera {
* <H2> Specification of Solution Thermodynamic Properties </H2>
* <HR>
*
* 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
/*!

145
Cantera/src/thermo/MolalityVPSSTP.cpp
View file

@ -24,6 +24,7 @@
#include "MolalityVPSSTP.h"
#include <iomanip>
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<std::string> pNames;
vector<double*> 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) {
;
}
}
}

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

@ -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

4
Cantera/src/thermo/PDSS_SSVol.cpp
View file

@ -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());

139
Cantera/src/thermo/PureFluidPhase.cpp
View file

@ -17,6 +17,7 @@
#include "../../../ext/tpx/utils.h"
#include <cstdlib>
#include <iomanip>
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<double> x(kk, 0.0);
std::vector<double> y(kk, 0.0);
std::vector<double> mu(kk, 0.0);
std::vector<double> a(kk, 0.0);
std::vector<double> ac(kk, 0.0);
std::vector<double> hbar(kk, 0.0);
std::vector<double> sbar(kk, 0.0);
std::vector<double> ubar(kk, 0.0);
std::vector<double> cpbar(kk, 0.0);
std::vector<double> vbar(kk, 0.0);
vector<std::string> pNames;
vector<double*> 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

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

@ -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

4
Cantera/src/thermo/State.cpp
View file

@ -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;

3
Cantera/src/thermo/State.h
View file

@ -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.

156
Cantera/src/thermo/ThermoPhase.cpp
View file

@ -21,6 +21,7 @@
#endif
#include "ThermoPhase.h"
#include <iomanip>
//@{
#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<std::string> pNames;
std::vector<doublereal *> 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) {
;
}
}
}

54
Cantera/src/thermo/ThermoPhase.h
View file

@ -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:

4
Cantera/src/thermo/VPSSMgr.cpp
View file

@ -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);

32
Cantera/src/thermo/VPStandardStateTP.h
View file

@ -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");
}

1
Cantera/src/thermo/WaterProps.cpp
View file

@ -501,6 +501,7 @@ namespace Cantera {
0.978197,
0.579829,
-0.202354};
//! parameter
const doublereal Hij[6][7] =
{