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

[Input] Create HMWSoln objects from YAML definitions

Browse source
This commit is contained in:
Ray Speth 2019-01-18 00:45:05 -05:00
parent 7d964dae47
commit 56612115f3
4 changed files with 227 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

100
src/thermo/HMWSoln.cpp
View file

@ -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<double>()) {
// 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<double>(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<AnyMap>();
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<string>()
&& 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<AnyMap>()) {
auto& species = item["species"].asVector<string>(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<AnyMap>();
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;

1
src/thermo/ThermoFactory.cpp
View file

@ -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";

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

@ -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

30
test/thermo/thermoFromYaml.cpp
View file

@ -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<MolalityVPSSTP*>(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);
}
}