Added excess volume capability for GibbsExcessVPSSTP and children
Corrected activity coefficients calculation for >2 species
This commit is contained in:
parent
f9d90adbcf
commit
d9b4031614
12 changed files with 895 additions and 65 deletions
|
|
@ -65,7 +65,8 @@ 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;
|
||||
|
|
@ -156,7 +157,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 +168,7 @@ namespace Cantera {
|
|||
}
|
||||
doublereal dd = meanMolecularWeight() / vtotal;
|
||||
State::setDensity(dd);
|
||||
delete [] vbar;
|
||||
}
|
||||
|
||||
void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p) {
|
||||
|
|
@ -320,7 +324,8 @@ 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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -301,6 +301,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 +334,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");
|
||||
}
|
||||
|
||||
//@}
|
||||
|
|
@ -562,7 +601,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;
|
||||
|
|
|
|||
|
|
@ -197,7 +197,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 +333,13 @@ namespace Cantera {
|
|||
void IonsFromNeutralVPSSTP::getActivityConcentrations(doublereal* c) const {
|
||||
getActivities(c);
|
||||
}
|
||||
|
||||
void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs,vector_fp& charges){
|
||||
coeffs = fm_neutralMolec_ions_;
|
||||
charges = m_speciesCharge;
|
||||
//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 +606,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();
|
||||
|
|
@ -717,7 +750,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 +824,129 @@ 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;
|
||||
|
||||
//! 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
|
||||
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 +970,7 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
|
||||
|
||||
void IonsFromNeutralVPSSTP::setMassFractions(const doublereal* const y) {
|
||||
GibbsExcessVPSSTP::setMassFractions(y);
|
||||
calcNeutralMoleculeMoleFractions();
|
||||
|
|
@ -836,7 +992,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 +1186,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 ,
|
||||
|
|
@ -1243,7 +1400,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 +1430,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 +1529,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 +1562,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 +1570,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 +1587,77 @@ 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;
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/*
|
||||
* 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;
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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,65 @@ 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;
|
||||
|
||||
|
||||
virtual void getDissociationCoeffs(vector_fp& coeffs, vector_fp& charges);
|
||||
|
||||
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;}
|
||||
|
||||
|
||||
//@}
|
||||
|
|
@ -680,6 +741,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 +756,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:
|
||||
|
|
@ -818,8 +896,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_;
|
||||
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -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,6 +1055,38 @@ 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];
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -565,6 +565,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 +591,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
|
||||
|
|
@ -608,11 +633,12 @@ 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
|
||||
*/
|
||||
virtual void getdlnActCoeffdlnC(doublereal *dlnActCoeffdlnC) const;
|
||||
virtual void getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const;
|
||||
virtual void getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const;
|
||||
|
||||
|
||||
//@}
|
||||
|
|
@ -761,7 +787,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 +837,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
|
||||
/*!
|
||||
|
|
|
|||
|
|
@ -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());
|
||||
|
|
|
|||
|
|
@ -196,6 +196,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;
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -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");
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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");
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue