From 8b47274765e4cfe15f489bd738fad0ed41fee69f Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 16 Jan 2019 16:35:51 -0500 Subject: [PATCH] [Input] Create DebyeHuckel and PDSS_Water objects from YAML definitions --- include/cantera/thermo/DebyeHuckel.h | 4 +- src/thermo/DebyeHuckel.cpp | 79 ++++++++++++++++------ src/thermo/PDSSFactory.cpp | 1 + src/thermo/Species.cpp | 4 +- src/thermo/ThermoFactory.cpp | 1 + test/data/thermo-models.yaml | 99 ++++++++++++++++++++++++++++ test/thermo/phaseConstructors.cpp | 12 ++-- test/thermo/thermoFromYaml.cpp | 53 +++++++++++++++ 8 files changed, 224 insertions(+), 29 deletions(-) diff --git a/include/cantera/thermo/DebyeHuckel.h b/include/cantera/thermo/DebyeHuckel.h index f9391f686..77c5292a7 100644 --- a/include/cantera/thermo/DebyeHuckel.h +++ b/include/cantera/thermo/DebyeHuckel.h @@ -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 diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp index 406f33add..ef856c7c5 100644 --- a/src/thermo/DebyeHuckel.cpp +++ b/src/thermo/DebyeHuckel.cpp @@ -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(); + setDebyeHuckelModel(node["model"].asString()); + if (node.hasKey("A_Debye")) { + if (node["A_Debye"].is() + && 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()) { + auto& species = item["species"].asVector(2); + setBeta(species[0], species[1], item["beta"].asDouble()); + } + } + } + // Solvent m_waterSS = dynamic_cast(providePDSS(0)); if (m_waterSS) { @@ -697,11 +741,8 @@ bool DebyeHuckel::addSpecies(shared_ptr 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 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 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); } diff --git a/src/thermo/PDSSFactory.cpp b/src/thermo/PDSSFactory.cpp index 3a57c9fa0..2c7b30074 100644 --- a/src/thermo/PDSSFactory.cpp +++ b/src/thermo/PDSSFactory.cpp @@ -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(); }); diff --git a/src/thermo/Species.cpp b/src/thermo/Species.cpp index 45479ca63..dd621fc28 100644 --- a/src/thermo/Species.cpp +++ b/src/thermo/Species.cpp @@ -58,11 +58,11 @@ shared_ptr 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 diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index 1ce39302d..e1bdb1be1 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -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(); }); diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml index dbea283d5..0fd951d84 100644 --- a/test/data/thermo-models.yaml +++ b/test/data/thermo-models.yaml @@ -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 diff --git a/test/thermo/phaseConstructors.cpp b/test/thermo/phaseConstructors.cpp index 88b2fb6db..1b48b726a 100644 --- a/test/thermo/phaseConstructors.cpp +++ b/test/thermo/phaseConstructors.cpp @@ -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); } diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp index 07428629e..4973cbb11 100644 --- a/test/thermo/thermoFromYaml.cpp +++ b/test/thermo/thermoFromYaml.cpp @@ -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(*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(*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); + } +}