diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index cc602add3..4c09ccc16 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -657,10 +657,108 @@ void HMWSoln::setCroppingCoefficients(double ln_gamma_k_min, CROP_ln_gamma_o_max = ln_gamma_o_max; } +vector_fp getSizedVector(const AnyMap& item, const std::string& key, size_t nCoeffs) +{ + vector_fp v; + if (item[key].is()) { + // Allow a single value to be given directly, rather than as a list of + // one item + v.push_back(item[key].asDouble()); + } else { + v = item[key].asVector(1, nCoeffs); + } + if (v.size() == 1 && nCoeffs == 5) { + // Adapt constant-temperature data to be compatible with the "complex" + // temperature model + v.resize(5, 0.0); + } + return v; +} + void HMWSoln::initThermo() { MolalityVPSSTP::initThermo(); - initLengths(); + if (m_input.hasKey("activity-data")) { + auto& actData = m_input["activity-data"].as(); + setPitzerTempModel(actData["temperature-model"].asString()); + initLengths(); + size_t nCoeffs = 1; + if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + nCoeffs = 2; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + nCoeffs = 5; + } + if (actData.hasKey("A_Debye")) { + if (actData["A_Debye"].is() + && actData["A_Debye"].asString() == "variable") { + setA_Debye(-1); + } else { + setA_Debye(actData.convert("A_Debye", "kg^0.5/gmol^0.5")); + } + } + if (actData.hasKey("max-ionic-strength")) { + setMaxIonicStrength(actData["max-ionic-strength"].asDouble()); + } + if (actData.hasKey("interactions")) { + for (auto& item : actData["interactions"].asVector()) { + auto& species = item["species"].asVector(1, 3); + size_t nsp = species.size(); + double q0 = charge(speciesIndex(species[0])); + double q1 = (nsp > 1) ? charge(speciesIndex(species[1])) : 0; + double q2 = (nsp == 3) ? charge(speciesIndex(species[2])) : 0; + if (nsp == 2 && q0 * q1 < 0) { + // Two species with opposite charges - binary salt + vector_fp beta0 = getSizedVector(item, "beta0", nCoeffs); + vector_fp beta1 = getSizedVector(item, "beta1", nCoeffs); + vector_fp beta2 = getSizedVector(item, "beta2", nCoeffs); + vector_fp Cphi = getSizedVector(item, "Cphi", nCoeffs); + if (beta0.size() != beta1.size() || beta0.size() != beta2.size() + || beta0.size() != Cphi.size()) { + throw CanteraError("HMWSoln::initThermo", "Inconsistent" + " binary salt array sizes ({}, {}, {}, {})", + beta0.size(), beta1.size(), beta2.size(), Cphi.size()); + } + double alpha1 = item["alpha1"].asDouble(); + double alpha2 = item.getDouble("alpha2", 0.0); + setBinarySalt(species[0], species[1], beta0.size(), + beta0.data(), beta1.data(), beta2.data(), Cphi.data(), + alpha1, alpha2); + } else if (nsp == 2 && q0 * q1 > 0) { + // Two species with like charges - "theta" interaction + vector_fp theta = getSizedVector(item, "theta", nCoeffs); + setTheta(species[0], species[1], theta.size(), theta.data()); + } else if (nsp == 2 && q0 * q1 == 0) { + // Two species, including at least one neutral + vector_fp lambda = getSizedVector(item, "lambda", nCoeffs); + setLambda(species[0], species[1], lambda.size(), lambda.data()); + } else if (nsp == 3 && q0 * q1 * q2 != 0) { + // Three charged species - "psi" interaction + vector_fp psi = getSizedVector(item, "psi", nCoeffs); + setPsi(species[0], species[1], species[2], + psi.size(), psi.data()); + } else if (nsp == 3 && q0 * q1 * q2 == 0) { + // Three species, including one neutral + vector_fp zeta = getSizedVector(item, "zeta", nCoeffs); + setZeta(species[0], species[1], species[2], + zeta.size(), zeta.data()); + } else if (nsp == 1) { + // single species (should be neutral) + vector_fp mu = getSizedVector(item, "mu", nCoeffs); + setMunnn(species[0], mu.size(), mu.data()); + } + } + } + if (actData.hasKey("cropping-coefficients")) { + auto& crop = actData["cropping-coefficients"].as(); + setCroppingCoefficients( + crop.getDouble("ln_gamma_k_min", -5.0), + crop.getDouble("ln_gamma_k_max", 15.0), + crop.getDouble("ln_gamma_o_min", -6.0), + crop.getDouble("ln_gamma_o_max", 3.0)); + } + } else { + initLengths(); + } for (int i = 0; i < 17; i++) { elambda[i] = 0.0; diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index cd7df820c..7eaa7dd5f 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -61,6 +61,7 @@ ThermoFactory::ThermoFactory() reg("LatticeSolid", []() { return new LatticeSolidPhase(); }); reg("Lattice", []() { return new LatticePhase(); }); reg("HMW", []() { return new HMWSoln(); }); + m_synonyms["HMW-electrolyte"] = "HMW"; reg("IdealSolidSolution", []() { return new IdealSolidSolnPhase(); }); reg("DebyeHuckel", []() { return new DebyeHuckel(); }); m_synonyms["Debye-Huckel"] = "DebyeHuckel"; diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml index b97b0495c..503e89657 100644 --- a/test/data/thermo-models.yaml +++ b/test/data/thermo-models.yaml @@ -119,6 +119,41 @@ phases: state: {T: 298, P: 1 atm, X: {H(s): 0.3, He(s): 0.7}} product-species: H(s) +- name: HMW-NaCl-electrolyte + species: [{HMW-species: [H2O(L), Cl-, H+, Na+, OH-]}] + thermo: HMW-electrolyte + state: + T: 423.15 + P: 101325 + molalities: {Na+: 6.0997, Cl-: 6.0996986044628, H+: 2.1628e-9, OH-: 1.3977e-6} + activity-data: + temperature-model: complex # "constant" or "linear" are the other options + A_Debye: 1.175930 kg^0.5/gmol^0.5 + interactions: + - species: [Na+, Cl-] + beta0: [0.0765, 0.008946, -3.3158E-6, -777.03, -4.4706] + beta1: [0.2664, 6.1608E-5, 1.0715E-6, 0.0, 0.0] + beta2: [0.0, 0.0, 0.0, 0.0, 0.0] + Cphi: [0.00127, -4.655E-5, 0.0, 33.317, 0.09421] + alpha1: 2.0 + - species: [H+, Cl-] + beta0: [0.1775] + beta1: [0.2945] + beta2: [0.0] + Cphi: [0.0008] + alpha1: 2.0 + - species: [Na+, OH-] + beta0: 0.0864 + beta1: 0.253 + beta2: 0.0 + Cphi: 0.0044 + alpha1: 2.0 + alpha2: 0.0 + - {species: [Cl-, OH-], theta: -0.05} + - {species: [Na+, Cl-, OH-], psi: -0.006} + - {species: [Na+, H+], theta: 0.036} + - {species: [Cl-, Na+, H+], psi: [-0.004]} + species: - name: NaCl(s) @@ -378,3 +413,65 @@ gas-species: -1063.94356, 3.65767573] - [3.28253784, 0.00148308754, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, -1088.45772, 5.45323129] + + +HMW-species: +- name: H2O(L) + composition: {H: 2, O: 1} + thermo: + model: NASA7 + temperature-ranges: [273.15, 600] + data: + - [7.255750050E+01, -6.624454020E-01, 2.561987460E-03, -4.365919230E-06, + 2.781789810E-09, -4.188654990E+04, -2.882801370E+02] + equation-of-state: + model: water-IAPWS95 + +- name: Na+ + composition: {Na: 1, E: -1} + thermo: + model: piecewise-Gibbs + reference-pressure: 1 bar + h0: 0.0 + dimensionless: true + data: {298.15: -125.5213, 333.15: -125.5213} + equation-of-state: + model: constant-volume + molar-volume: 1.3 + +- name: Cl- + composition: {Cl: 1, E: 1} + thermo: + model: piecewise-Gibbs + reference-pressure: 1 bar + h0: 0.0 + dimensionless: true + data: {298.15: -52.8716, 333.15: -52.8716} + equation-of-state: + reference-pressure: 1 bar + model: constant-volume + molar-volume: 1.3 + +- name: H+ + composition: {H: 1, E: -1} + thermo: + model: piecewise-Gibbs + reference-pressure: 1 bar + h0: 0.0 + dimensionless: true + data: {298.15: 0, 333.15: 0} + equation-of-state: + model: constant-volume + molar-volume: 1.3 + +- name: OH- + composition: {O: 1, H: 1, E: 1} + thermo: + model: piecewise-Gibbs + reference-pressure: 1 bar + h0: 0.0 + dimensionless: true + data: {298.15: -91.523, 333.15: -91.523} + equation-of-state: + model: constant-volume + molar-volume: 1.3 diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp index 879fdaad0..5e69af887 100644 --- a/test/thermo/thermoFromYaml.cpp +++ b/test/thermo/thermoFromYaml.cpp @@ -272,3 +272,33 @@ TEST(ThermoFromYaml, MaskellSolidSoln) EXPECT_NEAR(chemPotentials[0], -4.989677478024063e6, 1e-6); EXPECT_NEAR(chemPotentials[1], 4.989677478024063e6 + 1000, 1e-6); } + +TEST(ThermoFromYaml, HMWSoln) +{ + AnyMap infile = AnyMap::fromYamlFile("thermo-models.yaml"); + auto phaseNodes = infile["phases"].asMap("name"); + auto thermo = newPhase(*phaseNodes.at("HMW-NaCl-electrolyte"), infile); + + size_t N = thermo->nSpecies(); + auto HMW = dynamic_cast(thermo.get()); + vector_fp acMol(N), mf(N), activities(N), moll(N), mu0(N); + thermo->getMoleFractions(mf.data()); + HMW->getMolalities(moll.data()); + HMW->getMolalityActivityCoefficients(acMol.data()); + thermo->getActivities(activities.data()); + thermo->getStandardChemPotentials(mu0.data()); + + double acMolRef[] = {0.9341, 1.0191, 3.9637, 1.0191, 0.4660}; + double mfRef[] = {0.8198, 0.0901, 0.0000, 0.0901, 0.0000}; + double activitiesRef[] = {0.7658, 6.2164, 0.0000, 6.2164, 0.0000}; + double mollRef[] = {55.5084, 6.0997, 0.0000, 6.0997, 0.0000}; + double mu0Ref[] = {-317.175788, -186.014558, 0.0017225, -441.615429, -322.000412}; // kJ/gmol + + for (size_t k = 0 ; k < N; k++) { + EXPECT_NEAR(acMol[k], acMolRef[k], 2e-4); + EXPECT_NEAR(mf[k], mfRef[k], 2e-4); + EXPECT_NEAR(activities[k], activitiesRef[k], 2e-4); + EXPECT_NEAR(moll[k], mollRef[k], 2e-4); + EXPECT_NEAR(mu0[k]/1e6, mu0Ref[k], 2e-6); + } +}