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

[Thermo] Remove unused data from HMWSoln

Browse source 
The model implemented does not use the weak acid / "electrolyte species type"
concept, so there is no reason to read data about this from the input file.
This commit is contained in:
Ray Speth 2017-08-05 12:57:01 -04:00
parent cfc3b728f5
commit 6ca782030d
3 changed files with 0 additions and 257 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

126
include/cantera/thermo/HMWSoln.h
View file

@ -168,93 +168,6 @@ class WaterProps;
* molality has defined units of gmol kg-1, and therefore the ionic strength has
* units of sqrt(gmol/kg).
*
* In some instances, from some authors, a different formulation is used for the
* ionic strength in the equations below. The different formulation is due to
* the possibility of the existence of weak acids and how association wrt to the
* weak acid equilibrium relation affects the calculation of the activity
* coefficients via the assumed value of the ionic strength.
*
* If we are to assume that the association reaction doesn't have an effect on
* the ionic strength, then we will want to consider the associated weak acid as
* in effect being fully dissociated, when we calculate an effective value for
* the ionic strength. We will call this calculated value, the stoichiometric
* ionic strength, \f$ I_s \f$, putting a subscript s to denote it from the more
* straightforward calculation of \f$ I \f$.
*
* \f[
* I_s = \frac{1}{2} \sum_k{m_k^s z_k^2}
* \f]
*
* Here, \f$ m_k^s \f$ is the value of the molalities calculated assuming that
* all weak acid-base pairs are in their fully dissociated states. This
* calculation may be simplified by considering that the weakly associated acid
* may be made up of two charged species, k1 and k2, each with their own
* charges, obeying the following relationship:
*
* \f[
* z_k = z_{k1} + z_{k2}
* \f]
* Then, we may only need to specify one charge value, say, \f$ z_{k1}\f$, the
* cation charge number, in order to get both numbers, since we have already
* specified \f$ z_k \f$ in the definition of original species. Then, the
* stoichiometric ionic strength may be calculated via the following formula.
*
* \f[
* I_s = \frac{1}{2} \left(\sum_{k,ions}{m_k z_k^2}+
* \sum_{k,weak_assoc}(m_k z_{k1}^2 + m_k z_{k2}^2) \right)
* \f]
*
* The specification of which species are weakly associated acids is made in the
* input file via the `stoichIsMods` XML block, where the charge for k1 is also
* specified. An example is given below:
*
* @code
* <stoichIsMods>
* NaCl(aq):-1.0
* </stoichIsMods>
* @endcode
*
* Because we need the concept of a weakly associated acid in order to calculated
* \f$ I_s \f$ we need to catalog all species in the phase. This is done using
* the following categories:
*
* - `cEST_solvent` Solvent species (neutral)
* - `cEST_chargedSpecies` Charged species (charged)
* - `cEST_weakAcidAssociated` Species which can break apart into charged species.
* It may or may not be charged. These may or
* may not be be included in the
* species solution vector.
* - `cEST_strongAcidAssociated` Species which always breaks apart into charged species.
* It may or may not be charged. Normally, these
* aren't included in the speciation vector.
* - `cEST_polarNeutral` Polar neutral species
* - `cEST_nonpolarNeutral` Non polar neutral species
*
* Polar and non-polar neutral species are differentiated, because some
* additions to the activity coefficient expressions distinguish between these
* two types of solutes. This is the so-called salt-out effect.
*
* The type of species is specified in the `electrolyteSpeciesType` XML block.
* Note, this is not considered a part of the specification of the standard
* state for the species, at this time. Therefore, this information is put under
* the `activityCoefficient` XML block. An example is given below
*
* @code
* <electrolyteSpeciesType>
* H2L(L):solvent
* H+:chargedSpecies
* NaOH(aq):weakAcidAssociated
* NaCl(aq):strongAcidAssociated
* NH3(aq):polarNeutral
* O2(aq):nonpolarNeutral
* </electrolyteSpeciesType>
* @endcode
*
* Much of the species electrolyte type information is inferred from other
* information in the input file. For example, as species which is charged is
* given the "chargedSpecies" default category. A neutral solute species is put
* into the "nonpolarNeutral" category by default.
*
* ### Specification of the Excess Gibbs Free Energy
*
* Pitzer's formulation may best be represented as a specification of the excess
@ -1785,18 +1698,6 @@ private:
*/
int m_formPitzerTemp;
//! Vector containing the electrolyte species type
/*!
* The possible types are:
* - solvent
* - Charged Species
* - weakAcidAssociated
* - strongAcidAssociated
* - polarNeutral
* - nonpolarNeutral
*/
vector_int m_electrolyteSpeciesType;
//! a_k = Size of the ionic species in the DH formulation. units = meters
vector_fp m_Aionic;
@ -1812,11 +1713,6 @@ private:
//! Reference Temperature for the Pitzer formulations.
double m_TempPitzerRef;
//! Stoichiometric ionic strength on the molality scale. This differs from
//! m_IionicMolality in the sense that associated salts are treated as
//! unassociated salts, when calculating the Ionic strength by this method.
mutable double m_IionicMolalityStoich;
public:
/**
* Form of the constant outside the Debye-Huckel term called A. It's
@ -1884,21 +1780,6 @@ private:
//! vector of size m_kk, used as a temporary holding area.
mutable vector_fp m_tmpV;
/**
* Stoichiometric species charge -> This is for calculations of the ionic
* strength which ignore ion-ion pairing into neutral molecules. The
* Stoichiometric species charge is the charge of one of the ion that would
* occur if the species broke into two charged ion pairs.
*
* NaCl -> m_speciesCharge_Stoich = -1;
* HSO4- -> H+ + SO42- = -2
* -> The other charge is calculated.
*
* For species that aren't ion pairs, its equal to the m_speciesCharge[]
* value.
*/
vector_fp m_speciesCharge_Stoich;
/**
* Array of 2D data used in the Pitzer/HMW formulation. Beta0_ij[i][j] is
* the value of the Beta0 coefficient for the ij salt. It will be nonzero
@ -2658,13 +2539,6 @@ private:
//! Calculate molality cut-off parameters
void calcMCCutoffParams_();
//! Utility function to assign an integer value from a string for the
//! ElectrolyteSpeciesType field.
/*!
* @param estString string name of the electrolyte species type
*/
static int interp_est(const std::string& estString);
public:
//! Turn on copious debug printing when this is true
mutable int m_debugCalc;

20
src/thermo/HMWSoln.cpp
View file

@ -29,7 +29,6 @@ HMWSoln::HMWSoln() :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_WATER),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
@ -76,7 +75,6 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_WATER),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
@ -120,7 +118,6 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_WATER),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
@ -611,10 +608,7 @@ double HMWSoln::AionicRadius(int k) const
void HMWSoln::initLengths()
{
// Resize lengths equal to the number of species in the phase.
m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
m_speciesSize.resize(m_kk);
m_speciesCharge_Stoich.resize(m_kk, 0.0);
m_Aionic.resize(m_kk, 0.0);
m_tmpV.resize(m_kk, 0.0);
m_molalitiesCropped.resize(m_kk, 0.0);
@ -743,20 +737,6 @@ void HMWSoln::s_update_lnMolalityActCoeff() const
// coefficient calculations.
calcMolalitiesCropped();
// Calculate the stoichiometric ionic charge. This isn't used in the Pitzer
// formulation.
m_IionicMolalityStoich = 0.0;
for (size_t k = 0; k < m_kk; k++) {
double z_k = charge(k);
double zs_k1 = m_speciesCharge_Stoich[k];
if (z_k == zs_k1) {
m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
} else {
double zs_k2 = z_k - zs_k1;
m_IionicMolalityStoich += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
}
}
// Update the temperature dependence of the pitzer coefficients and their
// derivatives
s_updatePitzer_CoeffWRTemp();

111
src/thermo/HMWSoln_input.cpp
View file

@ -24,27 +24,6 @@ using namespace std;
namespace Cantera
{
int HMWSoln::interp_est(const std::string& estString)
{
if (ba::iequals(estString, "solvent")) {
return cEST_solvent;
} else if (ba::iequals(estString, "chargedspecies")) {
return cEST_chargedSpecies;
} else if (ba::iequals(estString, "weakacidassociated")) {
return cEST_weakAcidAssociated;
} else if (ba::iequals(estString, "strongacidassociated")) {
return cEST_strongAcidAssociated;
} else if (ba::iequals(estString, "polarneutral")) {
return cEST_polarNeutral;
} else if (ba::iequals(estString, "nonpolarneutral")) {
return cEST_nonpolarNeutral;
}
int retn, rval;
if ((retn = sscanf(estString.c_str(), "%d", &rval)) != 1) {
return -1;
}
return rval;
}
void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
{
@ -1032,13 +1011,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
// water calculator.
m_waterProps.reset(new WaterProps(&dynamic_cast<PDSS_Water&>(*m_waterSS)));
// Fill in parameters for the calculation of the stoichiometric Ionic
// Strength. The default is that stoich charge is the same as the regular
// charge.
for (size_t k = 0; k < m_kk; k++) {
m_speciesCharge_Stoich[k] = charge(k);
}
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
XML_Node* acNodePtr = 0;
@ -1084,40 +1056,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
// First look at the species database. Look for the subelement
// "stoichIsMods" in each of the species SS databases.
std::vector<const XML_Node*> xspecies = speciesData();
for (size_t k = 0; k < m_kk; k++) {
size_t jmap = npos;
string kname = speciesName(k);
for (size_t j = 0; j < xspecies.size(); j++) {
const XML_Node& sp = *xspecies[j];
string jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
}
}
if (jmap != npos) {
const XML_Node& sp = *xspecies[jmap];
getOptionalFloat(sp, "stoichIsMods", m_speciesCharge_Stoich[k]);
}
}
// Now look at the activity coefficient database
if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) {
XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
map<string, string> msIs;
getMap(sIsNode, msIs);
for (const auto& b : msIs) {
size_t kk = speciesIndex(b.first);
if (kk != npos) {
double val = fpValue(b.second);
m_speciesCharge_Stoich[kk] = val;
}
}
}
// Loop through the children getting multiple instances of parameters
if (acNodePtr) {
for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
@ -1149,55 +1087,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
readXMLCroppingCoefficients(acNode);
}
// Fill in the vector specifying the electrolyte species type
//
// First fill in default values. Everything is either a charge species, a
// nonpolar neutral, or the solvent.
for (size_t k = 0; k < m_kk; k++) {
if (fabs(charge(k)) > 0.0001) {
m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
if (fabs(m_speciesCharge_Stoich[k] - charge(k)) > 0.0001) {
m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
}
} else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) {
m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
} else {
m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
}
}
m_electrolyteSpeciesType[0] = cEST_solvent;
// First look at the species database. Look for the subelement
// "stoichIsMods" in each of the species SS databases.
std::vector<const XML_Node*> xspecies = speciesData();
for (size_t k = 0; k < m_kk; k++) {
const XML_Node* spPtr = xspecies[k];
if (spPtr && spPtr->hasChild("electrolyteSpeciesType")) {
string est = getChildValue(*spPtr, "electrolyteSpeciesType");
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
throw CanteraError("HMWSoln::initThermoXML",
"Bad electrolyte type: " + est);
}
}
}
// Then look at the phase thermo specification
if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) {
XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
map<string, string> msEST;
getMap(ESTNode, msEST);
for (const auto& b : msEST) {
size_t kk = speciesIndex(b.first);
if (kk != npos) {
string est = b.second;
if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
throw CanteraError("HMWSoln::initThermoXML",
"Bad electrolyte type: " + est);
}
}
}
}
IMS_typeCutoff_ = 2;
if (IMS_typeCutoff_ == 2) {
calcIMSCutoffParams_();