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

[Input] Create DebyeHuckel and PDSS_Water objects from YAML definitions

Browse source
This commit is contained in:
Ray Speth 2019-01-16 16:35:51 -05:00
parent 07aadbce6c
commit 8b47274765
8 changed files with 224 additions and 29 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -913,8 +913,8 @@ 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'.
//! 'dilute-limit', 'B-dot-with-variable-a', 'B-dot-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

79
src/thermo/DebyeHuckel.cpp
View file

@ -307,15 +307,20 @@ static int interp_est(const std::string& estString)
{
if (caseInsensitiveEquals(estString, "solvent")) {
return cEST_solvent;
} else if (caseInsensitiveEquals(estString, "chargedspecies")) {
} else if (estString == "charged-species"
|| caseInsensitiveEquals(estString, "chargedspecies")) {
return cEST_chargedSpecies;
} else if (caseInsensitiveEquals(estString, "weakacidassociated")) {
} else if (estString == "weak-acid-associated"
|| caseInsensitiveEquals(estString, "weakacidassociated")) {
return cEST_weakAcidAssociated;
} else if (caseInsensitiveEquals(estString, "strongacidassociated")) {
} else if (estString == "strong-acid-associated"
|| caseInsensitiveEquals(estString, "strongacidassociated")) {
return cEST_strongAcidAssociated;
} else if (caseInsensitiveEquals(estString, "polarneutral")) {
} else if (estString == "polar-neutral"
|| caseInsensitiveEquals(estString, "polarneutral")) {
return cEST_polarNeutral;
} else if (caseInsensitiveEquals(estString, "nonpolarneutral")) {
} else if (estString == "nonpolar-neutral"
|| caseInsensitiveEquals(estString, "nonpolarneutral")) {
return cEST_nonpolarNeutral;
} else {
throw CanteraError("interp_est (DebyeHuckel)",
@ -324,16 +329,21 @@ static int interp_est(const std::string& estString)
}
void DebyeHuckel::setDebyeHuckelModel(const std::string& model) {
if (model == "" || caseInsensitiveEquals(model, "Dilute_limit")) {
if (model == ""
|| model == "dilute-limit"
|| caseInsensitiveEquals(model, "Dilute_limit")) {
m_formDH = DHFORM_DILUTE_LIMIT;
} else if (caseInsensitiveEquals(model, "Bdot_with_variable_a")) {
} else if (model == "B-dot-with-variable-a"
|| caseInsensitiveEquals(model, "Bdot_with_variable_a")) {
m_formDH = DHFORM_BDOT_AK;
} else if (caseInsensitiveEquals(model, "Bdot_with_common_a")) {
} else if (model == "B-dot-with-common-a"
|| caseInsensitiveEquals(model, "Bdot_with_common_a")) {
m_formDH = DHFORM_BDOT_ACOMMON;
} else if (caseInsensitiveEquals(model, "Beta_ij")) {
} else if (caseInsensitiveEquals(model, "beta_ij")) {
m_formDH = DHFORM_BETAIJ;
m_Beta_ij.resize(m_kk, m_kk, 0.0);
} else if (caseInsensitiveEquals(model, "Pitzer_with_Beta_ij")) {
} else if (model == "Pitzer-with-beta_ij"
|| caseInsensitiveEquals(model, "Pitzer_with_Beta_ij")) {
m_formDH = DHFORM_PITZER_BETAIJ;
m_Beta_ij.resize(m_kk, m_kk, 0.0);
} else {
@ -557,6 +567,40 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
void DebyeHuckel::initThermo()
{
MolalityVPSSTP::initThermo();
if (m_input.hasKey("activity-data")) {
auto& node = m_input["activity-data"].as<AnyMap>();
setDebyeHuckelModel(node["model"].asString());
if (node.hasKey("A_Debye")) {
if (node["A_Debye"].is<string>()
&& node["A_Debye"].asString() == "variable") {
setA_Debye(-1);
} else {
setA_Debye(node.convert("A_Debye", "kg^0.5/gmol^0.5"));
}
}
if (node.hasKey("B_Debye")) {
setB_Debye(node.convert("B_Debye", "kg^0.5/gmol^0.5/m"));
}
if (node.hasKey("max-ionic-strength")) {
setMaxIonicStrength(node["max-ionic-strength"].asDouble());
}
if (node.hasKey("use-Helgeson-fixed-form")) {
useHelgesonFixedForm(node["use-Helgeson-fixed-form"].asBool());
}
if (node.hasKey("default-ionic-radius")) {
setDefaultIonicRadius(node.convert("default-ionic-radius", "m"));
}
if (node.hasKey("B-dot")) {
setB_dot(node["B-dot"].asDouble());
}
if (node.hasKey("beta")) {
for (auto& item : node["beta"].asVector<AnyMap>()) {
auto& species = item["species"].asVector<string>(2);
setBeta(species[0], species[1], item["beta"].asDouble());
}
}
}
// Solvent
m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
if (m_waterSS) {
@ -697,11 +741,8 @@ bool DebyeHuckel::addSpecies(shared_ptr<Species> spec)
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
}
// NAN will be replaced with default value
m_Aionic.push_back(spec->input.convert("ionic-radius", "m", NAN));
// Guess electrolyte species type based on charge properties
int est = cEST_nonpolarNeutral;
@ -709,8 +750,8 @@ bool DebyeHuckel::addSpecies(shared_ptr<Species> spec)
if (fabs(spec->charge) > 0.0001) {
est = cEST_chargedSpecies;
}
if (spec->extra.hasKey("weak_acid_charge")) {
stoichCharge = spec->extra["weak_acid_charge"].asDouble();
if (spec->input.hasKey("weak-acid-charge")) {
stoichCharge = spec->input["weak-acid-charge"].asDouble();
if (fabs(stoichCharge - spec->charge) > 0.0001) {
est = cEST_weakAcidAssociated;
}
@ -722,8 +763,8 @@ bool DebyeHuckel::addSpecies(shared_ptr<Species> spec)
}
// Apply override of the electrolyte species type
if (spec->extra.hasKey("electrolyte_species_type")) {
est = interp_est(spec->extra["electrolyte_species_type"].asString());
if (spec->input.hasKey("electrolyte-species-type")) {
est = interp_est(spec->input["electrolyte-species-type"].asString());
}
m_electrolyteSpeciesType.push_back(est);
}

1
src/thermo/PDSSFactory.cpp
View file

@ -25,6 +25,7 @@ PDSSFactory::PDSSFactory()
m_synonyms["constant-volume"] = "constant-incompressible";
reg("water", []() { return new PDSS_Water(); });
m_synonyms["waterPDSS"] = m_synonyms["waterIAPWS"] = "water";
m_synonyms["water-IAPWS95"] = "water";
reg("ions-from-neutral", []() { return new PDSS_IonsFromNeutral(); });
m_synonyms["IonFromNeutral"] = "ions-from-neutral";
reg("temperature_polynomial", []() { return new PDSS_SSVol(); });

4
src/thermo/Species.cpp
View file

@ -58,11 +58,11 @@ shared_ptr<Species> newSpecies(const XML_Node& species_node)
// Extra data used for some electrolyte species
if (species_node.hasChild("stoichIsMods")) {
s->extra["weak_acid_charge"] = getFloat(species_node, "stoichIsMods");
s->input["weak-acid-charge"] = getFloat(species_node, "stoichIsMods");
}
if (species_node.hasChild("electrolyteSpeciesType")) {
s->extra["electrolyte_species_type"] = species_node.child("electrolyteSpeciesType").value();
s->input["electrolyte-species-type"] = species_node.child("electrolyteSpeciesType").value();
}
// Extra data optionally used by LatticePhase

1
src/thermo/ThermoFactory.cpp
View file

@ -63,6 +63,7 @@ ThermoFactory::ThermoFactory()
reg("HMW", []() { return new HMWSoln(); });
reg("IdealSolidSolution", []() { return new IdealSolidSolnPhase(); });
reg("DebyeHuckel", []() { return new DebyeHuckel(); });
m_synonyms["Debye-Huckel"] = "DebyeHuckel";
reg("IdealMolalSolution", []() { return new IdealMolalSoln(); });
m_synonyms["ideal-molal-solution"] = "IdealMolalSolution";
reg("IdealGasVPSS", []() { return new IdealSolnGasVPSS(); });

99
test/data/thermo-models.yaml
View file

@ -47,6 +47,42 @@ phases:
P: 1 atm
molalities: {CH4(aq): 0.01, H2S(aq): 0.03, CO2(aq): 0.1}
- name: debye-huckel-B-dot-ak
thermo: Debye-Huckel
species: [{dh-electrolyte-species: all}]
activity-data:
model: B-dot-with-variable-a
A_Debye: 1.172576 kg^0.5/gmol^0.5
B_Debye: 3.2864e9 kg^0.5/gmol^0.5/m
B-dot: 0.0410
max-ionic-strength: 50.0
default-ionic-radius: 4.0 angstrom
state:
T: 300 K
P: 1 atm
M: {Na+: 9.3549, Cl-: 9.3549, H+: 1.0499E-8, OH-: 1.3765E-6, NaCl(aq): 0.98492}
- name: debye-huckel-beta_ij
thermo: Debye-Huckel
species: [{dh-electrolyte-species: all}]
activity-data:
model: beta_ij
max-ionic-strength: 3.0
use-Helgeson-fixed-form: true
default-ionic-radius: 3.042843 angstrom
beta:
- species: [H+, Cl-]
beta: 0.27
- species: [Na+, Cl-]
beta: 0.15
- species: [Na+, OH-]
beta: 0.06
state:
T: 300 K
P: 1 atm
M: {Na+: 3.0, Cl-: 3.0, H+: 1.0499E-8, OH-: 1.3765E-6, NaCl(aq): 0.98492}
species:
- name: NaCl(s)
composition: {Na: 1, Cl: 1}
@ -148,3 +184,66 @@ ideal-molal-fake-species:
equation-of-state:
model: constant-volume
molar-volume: 0.1
dh-electrolyte-species:
- name: H2O(l)
composition: {H: 2, O: 1}
equation-of-state:
model: water-IAPWS95
- name: Na+
composition: {Na: 1, E: -1}
thermo:
model: piecewise-Gibbs
h0: -240.34 kJ/mol
dimensionless: true
data: {298.15: -103.98186, 333.15: -103.98186}
equation-of-state:
model: constant-volume
molar-volume: 1.3
ionic-radius: 4 Å
- name: Cl-
composition: {Cl: 1, E: 1}
thermo:
model: piecewise-Gibbs
h0: -167.08 kJ/mol
dimensionless: true
data: {298.15: -74.20664, 333.15: -74.20664}
equation-of-state:
model: constant-volume
molar-volume: 1.3
ionic-radius: 3.0 angstrom
- name: H+
composition: {H: 1, E: -1}
thermo:
model: piecewise-Gibbs
h0: 0
dimensionless: true
data: {298.15: 0, 333.15: 0}
equation-of-state:
model: constant-volume
molar-volume: 0.0
ionic-radius: 9e-10 m
- name: OH-
composition: {O: 1, H: 1, E: 1}
thermo:
model: piecewise-Gibbs
h0: -230.015 kJ/mol
dimensionless: true
data: {298.15: -91.50963, 333.15: -85.0}
equation-of-state:
model: constant-volume
molar-volume: 1.3
ionic-radius: 3.5 angstrom
- name: NaCl(aq)
composition: {Na: 1, Cl: 1}
thermo:
model: piecewise-Gibbs
h0: -96.03E3 cal/mol
dimensionless: true
data: {298.15: -174.5057463, 333.15: -174.5057463}
equation-of-state:
model: constant-volume
molar-volume: 1.3
electrolyte-species-type: weak-acid-associated
weak-acid-charge: -1.0

12
test/thermo/phaseConstructors.cpp
View file

@ -426,22 +426,22 @@ TEST(DebyeHuckel, fromScratch)
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;
sNa->input["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;
sCl->input["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;
sH->input["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;
sOH->input["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.0;
sNaCl->extra["electrolyte_species_type"] = "weakAcidAssociated";
sNaCl->input["weak-acid-charge"] = -1.0;
sNaCl->input["electrolyte-species-type"] = "weakAcidAssociated";
for (auto& s : {sH2O, sNa, sCl, sH, sOH, sNaCl}) {
p.addSpecies(s);
}

53
test/thermo/thermoFromYaml.cpp
View file

@ -1,6 +1,7 @@
#include "gtest/gtest.h"
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/thermo/Elements.h"
#include "cantera/thermo/MolalityVPSSTP.h"
using namespace Cantera;
@ -138,3 +139,55 @@ TEST(ThermoFromYaml, IdealMolalSoln)
EXPECT_NEAR(thermo->gibbs_mole(), -3.8986e7, 1e3);
EXPECT_NEAR(thermo->density(), 12.058, 1e-3);
}
TEST(ThermoFromYaml, DebyeHuckel_bdot_ak)
{
AnyMap infile = AnyMap::fromYamlFile("thermo-models.yaml");
auto phaseNodes = infile["phases"].asMap("name");
auto thermo = newPhase(*phaseNodes.at("debye-huckel-B-dot-ak"), infile);
// Regression test based on XML input file
EXPECT_EQ(thermo->type(), "DebyeHuckel");
EXPECT_NEAR(thermo->density(), 60.296, 1e-2);
EXPECT_NEAR(thermo->cp_mass(), 1.58213e5, 1e0);
EXPECT_NEAR(thermo->entropy_mass(), 4.04222e3, 1e-2);
vector_fp actcoeff(thermo->nSpecies());
vector_fp mu_ss(thermo->nSpecies());
auto& molphase = dynamic_cast<MolalityVPSSTP&>(*thermo);
molphase.getMolalityActivityCoefficients(actcoeff.data());
thermo->getStandardChemPotentials(mu_ss.data());
double act_ref[] = {0.849231, 1.18392, 0.990068, 1.69245, 1.09349, 1.0};
double mu_ss_ref[] = {-3.06816e+08, -2.57956e+08, -1.84117e+08, 0.0,
-2.26855e+08, -4.3292e+08};
for (size_t k = 0; k < thermo->nSpecies(); k++) {
EXPECT_NEAR(actcoeff[k], act_ref[k], 1e-5);
EXPECT_NEAR(mu_ss[k], mu_ss_ref[k], 1e3);
}
}
TEST(ThermoFromYaml, DebyeHuckel_beta_ij)
{
AnyMap infile = AnyMap::fromYamlFile("thermo-models.yaml");
auto phaseNodes = infile["phases"].asMap("name");
auto thermo = newPhase(*phaseNodes.at("debye-huckel-beta_ij"), infile);
// Regression test based on XML input file
EXPECT_EQ(thermo->type(), "DebyeHuckel");
EXPECT_NEAR(thermo->density(), 122.264, 1e-3);
EXPECT_NEAR(thermo->cp_mass(), 81262.8, 1e-1);
EXPECT_NEAR(thermo->entropy_mass(), 4022.27, 1e-2);
vector_fp actcoeff(thermo->nSpecies());
vector_fp mu_ss(thermo->nSpecies());
auto& molphase = dynamic_cast<MolalityVPSSTP&>(*thermo);
molphase.getMolalityActivityCoefficients(actcoeff.data());
thermo->getStandardChemPotentials(mu_ss.data());
double act_ref[] = {0.959912, 1.16955, 1.16955, 2.40275, 0.681552, 1.0};
double mu_ss_ref[] = {-3.06816e+08, -2.57956e+08, -1.84117e+08, 0,
-2.26855e+08, -4.3292e+08};
for (size_t k = 0; k < thermo->nSpecies(); k++) {
EXPECT_NEAR(actcoeff[k], act_ref[k], 1e-5);
EXPECT_NEAR(mu_ss[k], mu_ss_ref[k], 1e3);
}
}