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

[Thermo] Allow instantiation of DebyeHuckel without XML

Browse source
This commit is contained in:
Ray Speth 2017-03-02 00:22:27 -05:00
parent 974bbc7da4
commit bfdc2b9e1d
6 changed files with 262 additions and 197 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

21
include/cantera/thermo/DebyeHuckel.h
View file

@ -814,6 +814,7 @@ public:
*/
virtual bool addSpecies(shared_ptr<Species> spec);
virtual void initThermo();
virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
//! Return the Debye Huckel constant as a function of temperature
@ -916,11 +917,31 @@ public:
*/
double AionicRadius(int k = 0) const;
//! Set the DebyeHuckel parameterization form. Must be one of
//! 'dilute_limit', 'Bdot_with_variable_a', 'Bdot_with_common_a',
//! 'Beta_ij', or 'Pitzer_with_Beta_ij'.
void setDebyeHuckelModel(const std::string& form);
//! Returns the form of the Debye-Huckel parameterization used
int formDH() const {
return m_formDH;
}
//! Set the A_Debye parameter. If a negative value is provided, enables
//! calculation of A_Debye using the detailed water equation of state.
void setA_Debye(double A);
void setB_Debye(double B) { m_B_Debye = B; }
void setB_dot(double bdot);
void setMaxIonicStrength(double Imax) { m_maxIionicStrength = Imax; }
void useHelgesonFixedForm(bool mode=true) { m_useHelgesonFixedForm = mode; }
//! Set the default ionic radius [m] for each species
void setDefaultIonicRadius(double value);
//! Set the value for the beta interaction between species sp1 and sp2.
void setBeta(const std::string& sp1, const std::string& sp2, double value);
//! Returns a reference to M_Beta_ij
Array2D& get_Beta_ij() {
return m_Beta_ij;

4
include/cantera/thermo/Species.h
View file

@ -7,6 +7,7 @@
#define CT_SPECIES_H
#include "cantera/base/ct_defs.h"
#include "cantera/base/AnyMap.h"
namespace Cantera
{
@ -52,6 +53,9 @@ public:
//! Thermodynamic data for the species
shared_ptr<SpeciesThermoInterpType> thermo;
//! Extra data used for specific models
AnyMap extra;
};
//! Create a new Species object from a 'species' XML_Node.

354
src/thermo/DebyeHuckel.cpp
View file

@ -14,6 +14,7 @@
#include "cantera/thermo/DebyeHuckel.h"
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/thermo/PDSS_Water.h"
#include "cantera/thermo/PDSS_ConstVol.h"
#include "cantera/thermo/electrolytes.h"
#include "cantera/base/stringUtils.h"
#include "cantera/base/ctml.h"
@ -324,12 +325,80 @@ static int interp_est(const std::string& estString)
return cEST_polarNeutral;
} else if (ba::iequals(estString, "nonpolarneutral")) {
return cEST_nonpolarNeutral;
} else {
throw CanteraError("interp_est (DebyeHuckel)",
"Invalid electrolyte species type '{}'", estString);
}
int retn, rval;
if ((retn = sscanf(estString.c_str(), "%d", &rval)) != 1) {
return -1;
}
void DebyeHuckel::setDebyeHuckelModel(const std::string& model) {
if (model == "" || ba::iequals(model, "Dilute_limit")) {
m_formDH = DHFORM_DILUTE_LIMIT;
} else if (ba::iequals(model, "Bdot_with_variable_a")) {
m_formDH = DHFORM_BDOT_AK;
} else if (ba::iequals(model, "Bdot_with_common_a")) {
m_formDH = DHFORM_BDOT_ACOMMON;
} else if (ba::iequals(model, "Beta_ij")) {
m_formDH = DHFORM_BETAIJ;
m_Beta_ij.resize(m_kk, m_kk, 0.0);
} else if (ba::iequals(model, "Pitzer_with_Beta_ij")) {
m_formDH = DHFORM_PITZER_BETAIJ;
m_Beta_ij.resize(m_kk, m_kk, 0.0);
} else {
throw CanteraError("DebyeHuckel::setDebyeHuckelModel",
"Unknown model '{}'", model);
}
return rval;
}
void DebyeHuckel::setA_Debye(double A)
{
if (A < 0) {
m_form_A_Debye = A_DEBYE_WATER;
} else {
m_form_A_Debye = A_DEBYE_CONST;
m_A_Debye = A;
}
}
void DebyeHuckel::setB_dot(double bdot)
{
if (m_formDH == DHFORM_BETAIJ || m_formDH == DHFORM_DILUTE_LIMIT ||
m_formDH == DHFORM_PITZER_BETAIJ) {
throw CanteraError("DebyeHuckel::setB_dot",
"B_dot entry in the wrong DH form");
}
// Set B_dot parameters for charged species
for (size_t k = 0; k < nSpecies(); k++) {
if (fabs(charge(k)) > 0.0001) {
m_B_Dot[k] = bdot;
} else {
m_B_Dot[k] = 0.0;
}
}
}
void DebyeHuckel::setDefaultIonicRadius(double value)
{
for (size_t k = 0; k < m_kk; k++) {
if (std::isnan(m_Aionic[k])) {
m_Aionic[k] = value;
}
}
}
void DebyeHuckel::setBeta(const std::string& sp1, const std::string& sp2,
double value)
{
size_t k1 = speciesIndex(sp1);
if (k1 == npos) {
throw CanteraError("DebyeHuckel::setBeta", "Species '{}' not found", sp1);
}
size_t k2 = speciesIndex(sp2);
if (k2 == npos) {
throw CanteraError("DebyeHuckel::setBeta", "Species '{}' not found", sp2);
}
m_Beta_ij(k1, k2) = value;
m_Beta_ij(k2, k1) = value;
}
void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
@ -349,32 +418,14 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
XML_Node& thermoNode = phaseNode.child("thermo");
// Determine the form of the Debye-Huckel model, m_formDH. We will use this
// information to size arrays below.
// Determine the form of the Debye-Huckel model, m_formDH. We will use
// this information to size arrays below. If there is no XML node named
// "activityCoefficients", assume that we are doing the extreme dilute
// limit assumption
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& scNode = thermoNode.child("activityCoefficients");
m_formDH = DHFORM_DILUTE_LIMIT;
std::string formString = scNode.attrib("model");
if (formString != "") {
if (formString == "Dilute_limit") {
m_formDH = DHFORM_DILUTE_LIMIT;
} else if (formString == "Bdot_with_variable_a") {
m_formDH = DHFORM_BDOT_AK;
} else if (formString == "Bdot_with_common_a") {
m_formDH = DHFORM_BDOT_ACOMMON;
} else if (formString == "Beta_ij") {
m_formDH = DHFORM_BETAIJ;
} else if (formString == "Pitzer_with_Beta_ij") {
m_formDH = DHFORM_PITZER_BETAIJ;
} else {
throw CanteraError("DebyeHuckel::initThermoXML",
"Unknown standardConc model: " + formString);
}
}
setDebyeHuckelModel(thermoNode.child("activityCoefficients")["model"]);
} else {
// If there is no XML node named "activityCoefficients", assume
// that we are doing the extreme dilute limit assumption
m_formDH = DHFORM_DILUTE_LIMIT;
setDebyeHuckelModel("Dilute_limit");
}
// Reconcile the solvent name and index.
@ -411,68 +462,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
" should be first species");
}
// Now go get the specification of the standard states for species in the
// solution. This includes the molar volumes data blocks for incompressible
// species.
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB =
get_XML_NameID("speciesData", speciesList["datasrc"],
&phaseNode.root());
const vector<string>&sss = speciesNames();
for (size_t k = 0; k < m_kk; k++) {
XML_Node* s = speciesDB->findByAttr("name", sss[k]);
if (!s) {
throw CanteraError("DebyeHuckel::initThermoXML",
"Species Data Base " + sss[k] + " not found");
}
XML_Node* ss = s->findByName("standardState");
if (!ss) {
throw CanteraError("DebyeHuckel::initThermoXML",
"Species " + sss[k] +
" standardState XML block not found");
}
std::string modelString = ss->attrib("model");
if (modelString == "") {
throw CanteraError("DebyeHuckel::initThermoXML",
"Species " + sss[k] +
" standardState XML block model attribute not found");
}
if (k == 0) {
if (ba::iequals(modelString, "wateriapws") || ba::iequals(modelString, "real_water") ||
ba::iequals(modelString, "waterpdss")) {
// Initialize the water standard state model
m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
if (!m_waterSS) {
throw CanteraError("HMWSoln::installThermoXML",
"Dynamic cast to PDSS_Water failed");
}
// Fill in the molar volume of water (m3/kmol) at standard
// conditions to fill in the m_speciesSize entry with something
// reasonable.
m_waterSS->setState_TP(300., OneAtm);
double dens = m_waterSS->density();
double mw = m_waterSS->molecularWeight();
m_speciesSize[0] = mw / dens;
} else if (ba::iequals(modelString, "constant_incompressible")) {
m_speciesSize[k] = getFloat(*ss, "molarVolume", "toSi");
} else {
throw CanteraError("DebyeHuckel::initThermoXML",
"Solvent SS Model \"" + modelString +
"\" is not known");
}
} else {
if (!ba::iequals(modelString, "constant_incompressible")) {
throw CanteraError("DebyeHuckel::initThermoXML",
"Solute SS Model \"" + modelString +
"\" is not known");
}
m_speciesSize[k] = getFloat(*ss, "molarVolume", "toSI");
}
}
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
XML_Node* acNodePtr = 0;
@ -486,59 +475,34 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
string modelString = ss->attrib("model");
if (modelString != "") {
if (ba::iequals(modelString, "water")) {
m_form_A_Debye = A_DEBYE_WATER;
setA_Debye(-1);
} else {
throw CanteraError("DebyeHuckel::initThermoXML",
"A_Debye Model \"" + modelString +
"\" is not known");
}
} else {
m_A_Debye = getFloat(acNode, "A_Debye");
setA_Debye(getFloat(acNode, "A_Debye"));
}
}
// Initialize the water property calculator. It will share the internal
// eos water calculator.
if (m_form_A_Debye == A_DEBYE_WATER) {
m_waterProps.reset(new WaterProps(m_waterSS));
}
// Look for parameters for B_Debye
if (acNode.hasChild("B_Debye")) {
m_B_Debye = getFloat(acNode, "B_Debye");
setB_Debye(getFloat(acNode, "B_Debye"));
}
// Look for parameters for B_dot
if (acNode.hasChild("B_dot")) {
if (m_formDH == DHFORM_BETAIJ ||
m_formDH == DHFORM_DILUTE_LIMIT ||
m_formDH == DHFORM_PITZER_BETAIJ) {
throw CanteraError("DebyeHuckel:init",
"B_dot entry in the wrong DH form");
}
double bdot_common = getFloat(acNode, "B_dot");
// Set B_dot parameters for charged species
for (size_t k = 0; k < m_kk; k++) {
double z_k = charge(k);
if (fabs(z_k) > 0.0001) {
m_B_Dot[k] = bdot_common;
} else {
m_B_Dot[k] = 0.0;
}
}
setB_dot(getFloat(acNode, "B_dot"));
}
// Look for Parameters for the Maximum Ionic Strength
if (acNode.hasChild("maxIonicStrength")) {
m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
setMaxIonicStrength(getFloat(acNode, "maxIonicStrength"));
}
// Look for Helgeson Parameters
if (acNode.hasChild("UseHelgesonFixedForm")) {
m_useHelgesonFixedForm = true;
} else {
m_useHelgesonFixedForm = false;
}
useHelgesonFixedForm(acNode.hasChild("UseHelgesonFixedForm"));
// Look for parameters for the Ionic radius
if (acNode.hasChild("ionicRadius")) {
@ -551,11 +515,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
if (irNode.hasAttrib("default")) {
std::string ads = irNode.attrib("default");
double ad = fpValue(ads);
for (size_t k = 0; k < m_kk; k++) {
m_Aionic[k] = ad * Afactor;
}
setDefaultIonicRadius(Afactor * fpValue(irNode.attrib("default")));
}
// If the Debye-Huckel form is BDOT_AK, we can have separate values
@ -578,8 +538,10 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
// ignore the lack of agreement (HKM -> may be changed in the
// future).
for (const auto& b : m) {
size_t kk = speciesIndex(b.first);
m_Aionic[kk] = fpValue(b.second) * Afactor;
size_t k = speciesIndex(b.first);
if (k != npos) {
m_Aionic[k] = fpValue(b.second) * Afactor;
}
}
}
}
@ -590,7 +552,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
if (acNode.hasChild("DHBetaMatrix")) {
if (m_formDH == DHFORM_BETAIJ ||
m_formDH == DHFORM_PITZER_BETAIJ) {
m_Beta_ij.resize(m_kk, m_kk, 0.0);
XML_Node& irNode = acNode.child("DHBetaMatrix");
const vector<string>& sn = speciesNames();
getMatrixValues(irNode, sn, sn, m_Beta_ij, true, true);
@ -600,39 +561,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
// Fill in parameters for the calculation of the stoichiometric Ionic
// Strength. The default is that stoich charge is the same as the
// regular charge.
m_speciesCharge_Stoich.resize(m_kk, 0.0);
for (size_t k = 0; k < m_kk; k++) {
m_speciesCharge_Stoich[k] = m_speciesCharge[k];
}
// 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();
size_t jj = xspecies.size();
for (size_t k = 0; k < m_kk; k++) {
size_t jmap = npos;
std::string kname = speciesName(k);
for (size_t j = 0; j < jj; j++) {
const XML_Node& sp = *xspecies[j];
std::string jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
}
}
if (jmap != npos) {
const XML_Node& sp = *xspecies[jmap];
if (sp.hasChild("stoichIsMods")) {
double val = getFloat(sp, "stoichIsMods");
m_speciesCharge_Stoich[k] = val;
}
}
}
// Now look at the activity coefficient database
// Override stoichiometric Ionic Strength based on the phase definition
if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) {
XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
map<std::string, std::string> msIs;
@ -645,39 +574,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
// 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(m_speciesCharge[k]) > 0.0001) {
m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
if (fabs(m_speciesCharge_Stoich[k] - m_speciesCharge[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[m_indexSolvent] = 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++) {
std::string kname = speciesName(k);
const XML_Node* spPtr = xspecies[k];
if (spPtr && spPtr->hasChild("electrolyteSpeciesType")) {
std::string est = getChildValue(*spPtr, "electrolyteSpeciesType");
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
throw CanteraError("DebyeHuckel:initThermoXML",
"Bad electrolyte type: " + est);
}
}
}
// Then look at the phase thermo specification
// Override electrolyte species type based on the phase definition
if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) {
XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
map<std::string, std::string> msEST;
@ -699,6 +596,41 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
void DebyeHuckel::initThermo()
{
MolalityVPSSTP::initThermo();
// Solvent
m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
if (m_waterSS) {
m_waterSS->setState_TP(300., OneAtm);
double dens = m_waterSS->density();
double mw = m_waterSS->molecularWeight();
m_speciesSize[0] = mw / dens;
// Initialize the water property calculator. It will share the internal
// eos water calculator.
if (m_form_A_Debye == A_DEBYE_WATER) {
m_waterProps.reset(new WaterProps(m_waterSS));
}
} else if (dynamic_cast<PDSS_ConstVol*>(providePDSS(0))) {
m_speciesSize[0] = providePDSS(0)->molarVolume();
} else {
throw CanteraError("DebyeHuckel::initThermo", "Solvent standard state"
" model must be WaterIAPWS or constant_incompressible.");
}
// Solutes
for (size_t k = 1; k < nSpecies(); k++) {
PDSS_ConstVol* ss = dynamic_cast<PDSS_ConstVol*>(providePDSS(k));
if (ss) {
m_speciesSize[k] = ss->molarVolume();
} else {
throw CanteraError("DebyeHuckel::initThermo", "Solute standard"
" state model must be constant_incompressible.");
}
}
}
double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const
{
double T = temperature();
@ -810,15 +742,43 @@ bool DebyeHuckel::addSpecies(shared_ptr<Species> spec)
{
bool added = MolalityVPSSTP::addSpecies(spec);
if (added) {
m_electrolyteSpeciesType.push_back(cEST_polarNeutral);
m_speciesSize.push_back(0.0);
m_Aionic.push_back(0.0);
m_lnActCoeffMolal.push_back(0.0);
m_dlnActCoeffMolaldT.push_back(0.0);
m_d2lnActCoeffMolaldT2.push_back(0.0);
m_dlnActCoeffMolaldP.push_back(0.0);
m_B_Dot.push_back(0.0);
m_tmpV.push_back(0.0);
if (spec->extra.hasKey("ionic_radius")) {
m_Aionic.push_back(spec->extra["ionic_radius"].asDouble());
} else {
m_Aionic.push_back(NAN); // NAN will be replaced with default value
}
// Guess electrolyte species type based on charge properties
int est = cEST_nonpolarNeutral;
double stoichCharge = spec->charge;
if (fabs(spec->charge) > 0.0001) {
est = cEST_chargedSpecies;
}
if (spec->extra.hasKey("weak_acid_charge")) {
stoichCharge = spec->extra["weak_acid_charge"].asDouble();
if (fabs(stoichCharge - spec->charge) > 0.0001) {
est = cEST_weakAcidAssociated;
}
}
m_speciesCharge_Stoich.push_back(stoichCharge);
if (m_electrolyteSpeciesType.size() == 0) {
est = cEST_solvent; // species 0 is the solvent
}
// Apply override of the electrolyte species type
if (spec->extra.hasKey("electrolyte_species_type")) {
est = interp_est(spec->extra["electrolyte_species_type"].asString());
}
m_electrolyteSpeciesType.push_back(est);
}
return added;
}

9
src/thermo/Species.cpp
View file

@ -51,6 +51,15 @@ shared_ptr<Species> newSpecies(const XML_Node& species_node)
s->transport->validate(*s);
}
// Extra data used for some electrolyte species
if (species_node.hasChild("stoichIsMods")) {
s->extra["weak_acid_charge"] = getFloat(species_node, "stoichIsMods");
}
if (species_node.hasChild("electrolyteSpeciesType")) {
s->extra["electrolyte_species_type"] = species_node.child("electrolyteSpeciesType").value();
}
return s;
}

64
test/thermo/phaseConstructors.cpp
View file

@ -9,9 +9,11 @@
#include "cantera/thermo/IonsFromNeutralVPSSTP.h"
#include "cantera/thermo/IdealSolnGasVPSS.h"
#include "cantera/thermo/IdealMolalSoln.h"
#include "cantera/thermo/DebyeHuckel.h"
#include "cantera/thermo/NasaPoly2.h"
#include "cantera/thermo/ShomatePoly.h"
#include "cantera/thermo/IdealGasPhase.h"
#include "cantera/thermo/Mu0Poly.h"
#include "cantera/base/ctml.h"
#include "cantera/base/stringUtils.h"
#include <fstream>
@ -28,6 +30,16 @@ shared_ptr<Species> make_species(const std::string& name,
return species;
}
shared_ptr<Species> make_species(const std::string& name,
const std::string& composition, double h298,
double T1, double mu1, double T2, double mu2)
{
auto species = make_shared<Species>(name, parseCompString(composition));
double coeffs[] = {2, h298, T1, mu1*GasConstant*T1, T2, mu2*GasConstant*T2};
species->thermo.reset(new Mu0Poly(200, 3500, 101325, coeffs));
return species;
}
class FixedChemPotSstpConstructorTest : public testing::Test
{
};
@ -311,4 +323,56 @@ TEST(IdealMolalSoln, fromScratch)
EXPECT_NEAR(p.density(), 12.058, 1e-3);
}
TEST(DebyeHuckel, fromScratch)
{
DebyeHuckel p;
p.addUndefinedElements();
auto sH2O = make_species("H2O(l)", "H:2, O:1", h2oliq_nasa_coeffs);
auto sNa = make_species("Na+", "Na:1, E:-1", -240.34e6,
298.15, -103.98186, 333.15, -103.98186);
sNa->charge = 1;
sNa->extra["ionic_radius"] = 4.0e-10;
auto sCl = make_species("Cl-", "Cl:1, E:1", -167.08e6,
298.15, -74.20664, 333.15, -74.20664);
sCl->charge = -1;
sCl->extra["ionic_radius"] = 3.0e-10;
auto sH = make_species("H+", "H:1, E:-1", 0.0, 298.15, 0.0, 333.15, 0.0);
sH->charge = 1;
sH->extra["ionic_radius"] = 9.0e-10;
auto sOH = make_species("OH-", "O:1, H:1, E:1", -230.015e6,
298.15, -91.50963, 333.15, -85);
sOH->charge = -1;
sOH->extra["ionic_radius"] = 3.5e-10;
auto sNaCl = make_species("NaCl(aq)", "Na:1, Cl:1", -96.03e6*4.184,
298.15, -174.5057463, 333.15, -174.5057463);
sNaCl->extra["weak_acid_charge"] = -1;
sNaCl->extra["electrolyte_species_type"] = "weakAcidAssociated";
for (auto& s : {sH2O, sNa, sCl, sH, sOH, sNaCl}) {
p.addSpecies(s);
}
size_t k = 0;
for (double v : {0.0555555, 0.0, 1.3, 1.3, 1.3, 1.3}) {
std::unique_ptr<PDSS_ConstVol> ss(new PDSS_ConstVol());
ss->setMolarVolume(v);
p.installPDSS(k++, std::move(ss));
}
p.setDebyeHuckelModel("bdot_with_variable_a");
p.setA_Debye(1.172576);
p.setB_Debye(3.2864e9);
p.setDefaultIonicRadius(3.5e-10);
p.setMaxIonicStrength(3.0);
p.useHelgesonFixedForm();
p.initThermo();
p.setState_TPM(300, 101325, "Na+:9.3549, Cl-:9.3549, H+:1.0499E-8,"
"OH-:1.3765E-6,NaCl(aq):0.98492");
// Regression test based on XML input file
vector_fp actcoeff(p.nSpecies());
p.getMolalityActivityCoefficients(actcoeff.data());
double act_ref[] = {1.21762, 0.538061, 0.472329, 0.717707, 0.507258, 1.0};
for (size_t k = 0; k < p.nSpecies(); k++) {
EXPECT_NEAR(actcoeff[k], act_ref[k], 1e-5);
}
}
} // namespace Cantera

7
test/thermo/thermo_data.h
View file

@ -6,6 +6,13 @@ const double h2o_nasa_coeffs[] = {
-9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+00,
4.19864056E+00, -2.03643410E-03, 6.52040211E-06, -5.48797062E-09,
1.77197817E-12, -3.02937267E+04, -8.49032208E-01};
const double h2oliq_nasa_coeffs[] = {
600.0, 7.255750050E+01, -6.624454020E-01, 2.561987460E-03,
-4.365919230E-06, 2.781789810E-09, -4.188671E+04, -2.8827879E+02,
7.255750050E+01, -6.624454020E-01, 2.561987460E-03, -4.365919230E-06,
2.781789810E-09, -4.188671E+04, -2.8827879E+02};
const double h2o_comp[] = {2.0, 1.0, 0.0};
const double h2_nasa_coeffs[] = {