diff --git a/include/cantera/thermo/LatticePhase.h b/include/cantera/thermo/LatticePhase.h index f7cbb09d8..69a0ed604 100644 --- a/include/cantera/thermo/LatticePhase.h +++ b/include/cantera/thermo/LatticePhase.h @@ -606,6 +606,8 @@ public: //! Set the density of lattice sites [kmol/m^3] void setSiteDensity(double sitedens); + virtual void initThermo(); + //! Set equation of state parameter values from XML entries. /*! * This method is called by function importPhase() when processing a phase diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp index c1a27862b..b73d0fe41 100644 --- a/src/thermo/LatticePhase.cpp +++ b/src/thermo/LatticePhase.cpp @@ -237,11 +237,26 @@ bool LatticePhase::addSpecies(shared_ptr spec) m_g0_RT.push_back(0.0); m_cp0_R.push_back(0.0); m_s0_R.push_back(0.0); - if (spec->extra.hasKey("molar_volume")) { - m_speciesMolarVolume.push_back(spec->extra["molar_volume"].asDouble()); - } else { - m_speciesMolarVolume.push_back(1.0 / m_site_density); + double mv = 1.0 / m_site_density; + if (spec->input.hasKey("equation-of-state")) { + auto& eos = spec->input["equation-of-state"].as(); + if (eos.getString("model", "") != "constant-volume") { + throw CanteraError("LatticePhase::initThermo", + "lattice model requires constant-volume species model " + "for species '{}'", spec->name); + } + if (eos.hasKey("density")) { + mv = molecularWeight(m_kk-1) / eos.convert("density", "kg/m^3"); + } else if (eos.hasKey("molar-density")) { + mv = 1.0 / eos.convert("molar-density", "kmol/m^3"); + } else if (eos.hasKey("molar-volume")) { + mv = eos.convert("molar-volume", "m^3/kmol"); + } + } else if (spec->extra.hasKey("molar_volume")) { + // from XML + mv = spec->extra["molar_volume"].asDouble(); } + m_speciesMolarVolume.push_back(mv); } return added; } @@ -249,6 +264,18 @@ bool LatticePhase::addSpecies(shared_ptr spec) void LatticePhase::setSiteDensity(double sitedens) { m_site_density = sitedens; + for (size_t k = 0; k < m_kk; k++) { + if (species(k)->extra.hasKey("molar_volume")) { + continue; + } else if (species(k)->input.hasKey("equation-of-state")) { + auto& eos = species(k)->input["equation-of-state"]; + if (eos.hasKey("molar-volume") || eos.hasKey("density") + || eos.hasKey("molar-density")) { + continue; + } + } + m_speciesMolarVolume[k] = 1.0 / m_site_density; + } } void LatticePhase::_updateThermo() const @@ -264,6 +291,13 @@ void LatticePhase::_updateThermo() const } } +void LatticePhase::initThermo() +{ + if (m_input.hasKey("site-density")) { + setSiteDensity(m_input.convert("site-density", "kmol/m^3")); + } +} + void LatticePhase::setParametersFromXML(const XML_Node& eosdata) { eosdata._require("model", "Lattice"); diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp index f7451a338..771ad12a0 100644 --- a/src/thermo/LatticeSolidPhase.cpp +++ b/src/thermo/LatticeSolidPhase.cpp @@ -288,46 +288,19 @@ void LatticeSolidPhase::getGibbs_ref(doublereal* g) const void LatticeSolidPhase::initThermo() { - size_t kk = 0; - size_t kstart = 0; - lkstart_.resize(m_lattice.size() + 1); - size_t loc = 0; - - for (size_t n = 0; n < m_lattice.size(); n++) { - shared_ptr& lp = m_lattice[n]; - vector_fp constArr(lp->nElements()); - const vector_fp& aws = lp->atomicWeights(); - for (size_t es = 0; es < lp->nElements(); es++) { - addElement(lp->elementName(es), aws[es], lp->atomicNumber(es), - lp->entropyElement298(es), lp->elementType(es)); - } - kstart = kk; - - for (size_t k = 0; k < lp->nSpecies(); k++) { - addSpecies(lp->species(k)); - kk++; - } - // Add in the lattice stoichiometry constraint - if (n > 0) { - string econ = fmt::format("LC_{}_{}", n, id()); - size_t m = addElement(econ, 0.0, 0, 0.0, CT_ELEM_TYPE_LATTICERATIO); - size_t mm = nElements(); - size_t nsp0 = m_lattice[0]->nSpecies(); - for (size_t k = 0; k < nsp0; k++) { - m_speciesComp[k * mm + m] = -theta_[0]; - } - for (size_t k = 0; k < lp->nSpecies(); k++) { - size_t ks = kstart + k; - m_speciesComp[ks * mm + m] = theta_[n]; + if (m_input.hasKey("composition") && m_input.hasKey("__file__")) { + AnyMap infile = AnyMap::fromYamlFile(m_input["__file__"].asString()); + auto phaseNodes = infile["phases"].asMap("name"); + compositionMap composition = m_input["composition"].asMap(); + for (auto& item : composition) { + if (phaseNodes.count(item.first)) { + addLattice(newPhase(*phaseNodes.at(item.first), infile)); + } else { + throw CanteraError("LatticeSolidPhase::initThermo", + "Could not find component phase named '{}'.", item.first); } } - size_t nsp = m_lattice[n]->nSpecies(); - lkstart_[n] = loc; - for (size_t k = 0; k < nsp; k++) { - m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_lattice.size(); - loc++; - } - lkstart_[n+1] = loc; + setLatticeStoichiometry(composition); } setMoleFractions(m_x.data()); @@ -336,22 +309,35 @@ void LatticeSolidPhase::initThermo() bool LatticeSolidPhase::addSpecies(shared_ptr spec) { - bool added = ThermoPhase::addSpecies(spec); - if (added) { - m_x.push_back(0.0); - tmpV_.push_back(0.0); - } - return added; + // Species are added from component phases in addLattice() + return false; } void LatticeSolidPhase::addLattice(shared_ptr lattice) { m_lattice.push_back(lattice); + if (lkstart_.empty()) { + lkstart_.push_back(0); + } + lkstart_.push_back(lkstart_.back() + lattice->nSpecies()); + if (theta_.size() == 0) { theta_.push_back(1.0); } else { theta_.push_back(0.0); } + + for (size_t k = 0; k < lattice->nSpecies(); k++) { + ThermoPhase::addSpecies(lattice->species(k)); + vector_fp constArr(lattice->nElements()); + const vector_fp& aws = lattice->atomicWeights(); + for (size_t es = 0; es < lattice->nElements(); es++) { + addElement(lattice->elementName(es), aws[es], lattice->atomicNumber(es), + lattice->entropyElement298(es), lattice->elementType(es)); + } + m_x.push_back(lattice->moleFraction(k)); + tmpV_.push_back(0.0); + } } void LatticeSolidPhase::setLatticeStoichiometry(const compositionMap& comp) @@ -359,6 +345,19 @@ void LatticeSolidPhase::setLatticeStoichiometry(const compositionMap& comp) for (size_t i = 0; i < m_lattice.size(); i++) { theta_[i] = getValue(comp, m_lattice[i]->name(), 0.0); } + // Add in the lattice stoichiometry constraint + for (size_t i = 1; i < m_lattice.size(); i++) { + string econ = fmt::format("LC_{}_{}", i, id()); + size_t m = addElement(econ, 0.0, 0, 0.0, CT_ELEM_TYPE_LATTICERATIO); + size_t mm = nElements(); + for (size_t k = 0; k < m_lattice[0]->nSpecies(); k++) { + m_speciesComp[k * mm + m] = -theta_[0]; + } + for (size_t k = 0; k < m_lattice[i]->nSpecies(); k++) { + size_t ks = lkstart_[i] + k; + m_speciesComp[ks * mm + m] = theta_[i]; + } + } } void LatticeSolidPhase::_updateThermo() const diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index 6ad6777ae..aeeb44b6c 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -63,7 +63,9 @@ ThermoFactory::ThermoFactory() reg("PureFluid", []() { return new PureFluidPhase(); }); m_synonyms["pure-fluid"] = "PureFluid"; reg("LatticeSolid", []() { return new LatticeSolidPhase(); }); + m_synonyms["compound-lattice"] = "LatticeSolid"; reg("Lattice", []() { return new LatticePhase(); }); + m_synonyms["lattice"] = "Lattice"; reg("HMW", []() { return new HMWSoln(); }); m_synonyms["HMW-electrolyte"] = "HMW"; reg("IdealSolidSolution", []() { return new IdealSolidSolnPhase(); }); diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml index f3076f73b..5a2ae0efc 100644 --- a/test/data/thermo-models.yaml +++ b/test/data/thermo-models.yaml @@ -191,6 +191,18 @@ phases: species: [{ISSP-species: all}] state: {T: 500, P: 2 bar, X: {sp1: 0.1, sp2: 0.89, sp3: 0.01}} +- name: Li7Si3(s) + species: [{lattice-species: [Li7Si3(s)]}] + thermo: fixed-stoichiometry +- name: Li7Si3-interstitial + species: [{lattice-species: [Li(i), V(i)]}] + thermo: lattice + site-density: 1.046344e-2 gmol/cm^3 + state: {T: 725 K, P: 10 atm, X: {Li(i): 0.01, V(i): 0.99}} +- name: Li7Si3_and_interstitials + thermo: compound-lattice + composition: {Li7Si3(s): 1.0, Li7Si3-interstitial: 1.0} + state: {T: 725 K, P: 10 atm} species: - name: NaCl(s) @@ -653,3 +665,32 @@ ISSP-species: equation-of-state: model: constant-volume molar-volume: 0.1 + + +lattice-species: +- name: Li7Si3(s) + composition: {Li: 7, Si: 3} + thermo: + model: Shomate + temperature-ranges: [250, 700, 2700] + data: + - [295.73961, -6.753295, -44.538551, 29.738846, -1.022387, -348.88919, 554.35647] + - [250.51429, 51.125155, -28.341244, 6.242135, 1.346861, -328.46578, 498.84106] + equation-of-state: + model: constant-volume + density: 1.39 g/cm^3 +- name: Li(i) + composition: {Li: 1} + thermo: + model: constant-cp + T0: 298.15 + s0: 20 J/mol/K + cp0: 20 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 0.2 +- name: V(i) + composition: {} + thermo: + model: constant-cp + h0: 89.8 J/mol diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp index 8f894eaa2..4ba1a710d 100644 --- a/test/thermo/thermoFromYaml.cpp +++ b/test/thermo/thermoFromYaml.cpp @@ -348,3 +348,22 @@ TEST(ThermoFromYaml, IdealSolidSolnPhase) EXPECT_NEAR(thermo->enthalpy_mass(), -15642803.3884617, 1e-4); EXPECT_NEAR(thermo->gibbs_mole(), -313642293.1654253, 1e-4); } + +TEST(ThermoFromYaml, Lattice) +{ + auto thermo = newThermo("thermo-models.yaml", "Li7Si3_and_interstitials"); + + // Regression test based on modified version of Li7Si3_ls.xml + EXPECT_NEAR(thermo->enthalpy_mass(), -2077821.9295456698, 1e-6); + double mu_ref[] = {-4.62717474e+08, -4.64248485e+07, 1.16370186e+05}; + double vol_ref[] = {0.09557086, 0.2, 0.09557086}; + vector_fp mu(thermo->nSpecies()); + vector_fp vol(thermo->nSpecies()); + thermo->getChemPotentials(mu.data()); + thermo->getPartialMolarVolumes(vol.data()); + + for (size_t k = 0; k < thermo->nSpecies(); k++) { + EXPECT_NEAR(mu[k], mu_ref[k], 1e-7*fabs(mu_ref[k])); + EXPECT_NEAR(vol[k], vol_ref[k], 1e-7); + } +}