From f69ef446007ee73ad8c5f527605196b514feca6e Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 17 Jun 2017 21:57:13 -0400 Subject: [PATCH] [Thermo] Allow instantiation of MargulesVPSSTP without XML --- include/cantera/thermo/MargulesVPSSTP.h | 24 +++-- src/thermo/MargulesVPSSTP.cpp | 134 +++++++++++------------- test/thermo/phaseConstructors.cpp | 34 ++++++ test/thermo/thermo_data.h | 6 ++ 4 files changed, 116 insertions(+), 82 deletions(-) diff --git a/include/cantera/thermo/MargulesVPSSTP.h b/include/cantera/thermo/MargulesVPSSTP.h index 4592ea340..8f293c29d 100644 --- a/include/cantera/thermo/MargulesVPSSTP.h +++ b/include/cantera/thermo/MargulesVPSSTP.h @@ -356,6 +356,23 @@ public: virtual void initThermo(); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); + //! Add a binary species interaction with the specified parameters + /*! + * @param speciesA name of the first species + * @param speciesB name of the second species + * @param h0 first excess enthalpy coefficient + * @param h1 second excess enthalpy coefficient + * @param s0 first excess entropy coefficient + * @param s1 second excess entropy coefficient + * @param vh0 first enthalpy coefficient for excess volume + * @param vh1 second enthalpy coefficient for excess volume + * @param vs0 first entropy coefficient for excess volume + * @param vs1 second entropy coefficient for excess volume + */ + void addBinaryInteraction(const std::string& speciesA, + const std::string& speciesB, double h0, double h1, double s0, double s1, + double vh0, double vh1, double vs0, double vs1); + //! @} //! @name Derivatives of Thermodynamic Variables needed for Applications //! @{ @@ -380,13 +397,6 @@ private: */ void readXMLBinarySpecies(XML_Node& xmlBinarySpecies); - //! Resize internal arrays within the object that depend upon the number of - //! binary Margules interaction terms - /*! - * @param num Number of binary Margules interaction terms - */ - void resizeNumInteractions(const size_t num); - //! Initialize lengths of local variables after all species have been //! identified. void initLengths(); diff --git a/src/thermo/MargulesVPSSTP.cpp b/src/thermo/MargulesVPSSTP.cpp index ff3ceccb5..3a0e5a9b0 100644 --- a/src/thermo/MargulesVPSSTP.cpp +++ b/src/thermo/MargulesVPSSTP.cpp @@ -254,6 +254,32 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_) GibbsExcessVPSSTP::initThermoXML(phaseNode, id_); } +void MargulesVPSSTP::addBinaryInteraction(const std::string& speciesA, + const std::string& speciesB, double h0, double h1, double s0, double s1, + double vh0, double vh1, double vs0, double vs1) +{ + size_t kA = speciesIndex(speciesA); + size_t kB = speciesIndex(speciesB); + // The interaction is silently ignored if either species is not defined in + // the current phase. + if (kA == npos || kB == npos) { + return; + } + m_pSpecies_A_ij.push_back(kA); + m_pSpecies_B_ij.push_back(kB); + + m_HE_b_ij.push_back(h0); + m_HE_c_ij.push_back(h1); + m_SE_b_ij.push_back(s0); + m_SE_c_ij.push_back(s1); + m_VHE_b_ij.push_back(vh0); + m_VHE_c_ij.push_back(vh1); + m_VSE_b_ij.push_back(vs0); + m_VSE_c_ij.push_back(vs1); + numBinaryInteractions_++; +} + + void MargulesVPSSTP::s_update_lnActCoeff() const { double T = temperature(); @@ -468,22 +494,6 @@ void MargulesVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnActCoeff } } -void MargulesVPSSTP::resizeNumInteractions(const size_t num) -{ - numBinaryInteractions_ = num; - m_HE_b_ij.resize(num, 0.0); - m_HE_c_ij.resize(num, 0.0); - m_SE_b_ij.resize(num, 0.0); - m_SE_c_ij.resize(num, 0.0); - m_VHE_b_ij.resize(num, 0.0); - m_VHE_c_ij.resize(num, 0.0); - m_VSE_b_ij.resize(num, 0.0); - m_VSE_c_ij.resize(num, 0.0); - - m_pSpecies_A_ij.resize(num, npos); - m_pSpecies_B_ij.resize(num, npos); -} - void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies) { string xname = xmLBinarySpecies.name(); @@ -491,7 +501,6 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies) throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "Incorrect name for processing this routine: " + xname); } - vector_fp vParams; string aName = xmLBinarySpecies.attrib("speciesA"); if (aName == "") { throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "no speciesA attrib"); @@ -501,35 +510,15 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies) throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "no speciesB attrib"); } - // Find the index of the species in the current phase. It's not an error to - // not find the species. What this means is that the A-B interaction - // referred to in this block will be ignored. - size_t aSpecies = speciesIndex(aName); - if (aSpecies == npos) { - return; - } - string aspName = speciesName(aSpecies); - - // @TODO Figure out what the original reason is for putting an error - // condition for charged species. Seems OK to me. - if (charge(aSpecies) != 0.0) { - throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", - "speciesA has a charge: {}", charge(aSpecies)); - } - size_t bSpecies = speciesIndex(bName); - if (bSpecies == npos) { - return; - } - string bspName = speciesName(bSpecies); - if (charge(bSpecies) != 0.0) { - throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", - "speciesB has a charge: {}", charge(bSpecies)); - } - - resizeNumInteractions(numBinaryInteractions_ + 1); - size_t iSpot = numBinaryInteractions_ - 1; - m_pSpecies_A_ij[iSpot] = aSpecies; - m_pSpecies_B_ij[iSpot] = bSpecies; + vector_fp vParams; + double h0 = 0.0; + double h1 = 0.0; + double s0 = 0.0; + double s1 = 0.0; + double vh0 = 0.0; + double vh1 = 0.0; + double vs0 = 0.0; + double vs1 = 0.0; for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) { XML_Node& xmlChild = xmLBinarySpecies.child(iChild); @@ -547,50 +536,45 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies) // Get the string containing all of the values getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy"); if (vParams.size() != 2) { - throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEnthalpy for " + aspName - + "::" + bspName, - "wrong number of params found. Need 2"); + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies" + "excessEnthalpy for {} : {}: wrong number of params found." + " Need 2", aName, bName); } - m_HE_b_ij[iSpot] = vParams[0]; - m_HE_c_ij[iSpot] = vParams[1]; - } - - if (nodeName == "excessentropy") { + h0 = vParams[0]; + h1 = vParams[1]; + } else if (nodeName == "excessentropy") { // Get the string containing all of the values getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy"); if (vParams.size() != 2) { - throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEntropy for " + aspName - + "::" + bspName, - "wrong number of params found. Need 2"); + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies" + "excessEntropy for {} : {}: wrong number of params found." + " Need 2", aName, bName); } - m_SE_b_ij[iSpot] = vParams[0]; - m_SE_c_ij[iSpot] = vParams[1]; - } - - if (nodeName == "excessvolume_enthalpy") { + s0 = vParams[0]; + s1 = vParams[1]; + } else if (nodeName == "excessvolume_enthalpy") { // Get the string containing all of the values getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy"); if (vParams.size() != 2) { - throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Enthalpy for " + aspName - + "::" + bspName, - "wrong number of params found. Need 2"); + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies" + "excessVolume_Enthalpy for {} : {}: wrong number of params" + " found. Need 2", aName, bName); } - m_VHE_b_ij[iSpot] = vParams[0]; - m_VHE_c_ij[iSpot] = vParams[1]; - } - - if (nodeName == "excessvolume_entropy") { + vh0 = vParams[0]; + vh1 = vParams[1]; + } else if (nodeName == "excessvolume_entropy") { // Get the string containing all of the values getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy"); if (vParams.size() != 2) { - throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Entropy for " + aspName - + "::" + bspName, - "wrong number of params found. Need 2"); + throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies" + "excessVolume_Entropy for {} : {}: wrong number of params" + " found. Need 2", aName, bName); } - m_VSE_b_ij[iSpot] = vParams[0]; - m_VSE_c_ij[iSpot] = vParams[1]; + vs0 = vParams[0]; + vs1 = vParams[1]; } } + addBinaryInteraction(aName, bName, h0, h1, s0, s1, vh0, vh1, vs0, vs1); } } diff --git a/test/thermo/phaseConstructors.cpp b/test/thermo/phaseConstructors.cpp index ee32be41b..974763438 100644 --- a/test/thermo/phaseConstructors.cpp +++ b/test/thermo/phaseConstructors.cpp @@ -10,6 +10,7 @@ #include "cantera/thermo/IdealSolnGasVPSS.h" #include "cantera/thermo/IdealMolalSoln.h" #include "cantera/thermo/DebyeHuckel.h" +#include "cantera/thermo/MargulesVPSSTP.h" #include "cantera/thermo/NasaPoly2.h" #include "cantera/thermo/ShomatePoly.h" #include "cantera/thermo/IdealGasPhase.h" @@ -30,6 +31,14 @@ shared_ptr make_species(const std::string& name, return species; } +shared_ptr make_shomate_species(const std::string& name, + const std::string& composition, const double* shomate_coeffs) +{ + auto species = make_shared(name, parseCompString(composition)); + species->thermo.reset(new ShomatePoly(200, 3500, 101325, shomate_coeffs)); + return species; +} + shared_ptr make_species(const std::string& name, const std::string& composition, double h298, double T1, double mu1, double T2, double mu2) @@ -375,4 +384,29 @@ TEST(DebyeHuckel, fromScratch) } } +TEST(MargulesVPSSTP, fromScratch) +{ + MargulesVPSSTP p; + p.addUndefinedElements(); + auto sKCl = make_shomate_species("KCl(L)", "K:1 Cl:1", kcl_shomate_coeffs); + auto sLiCl = make_shomate_species("LiCl(L)", "Li:1 Cl:1", licl_shomate_coeffs); + p.addSpecies(sKCl); + p.addSpecies(sLiCl); + size_t k = 0; + for (double v : {0.03757, 0.020304}) { + std::unique_ptr ss(new PDSS_ConstVol()); + ss->setMolarVolume(v); + p.installPDSS(k++, std::move(ss)); + } + p.initThermo(); + p.setState_TPX(900, 101325, "KCl(L):0.3, LiCl(L):0.7"); + p.addBinaryInteraction("KCl(L)", "LiCl(L)", + -1.757e7, -3.77e5, -7.627e3, 4.958e3, 0.0, 0.0, 0.0, 0.0); + + // Regression test based on LiKCl_liquid.xml + EXPECT_NEAR(p.density(), 2042.1165603245981, 1e-9); + EXPECT_NEAR(p.gibbs_mass(), -9682981.421693124, 1e-5); + EXPECT_NEAR(p.cp_mole(), 67478.48085733457, 1e-8); +} + } // namespace Cantera diff --git a/test/thermo/thermo_data.h b/test/thermo/thermo_data.h index 36cb53148..21f63377d 100644 --- a/test/thermo/thermo_data.h +++ b/test/thermo/thermo_data.h @@ -47,4 +47,10 @@ const double co_shomate_coeffs[] = { 35.15070, 1.300095, -0.205921, 0.013550, -3.282780, -127.8375, 231.7120}; const double co_comp[] = {0.0, 1.0, 1.0}; +// single-region Shomate coefficients +const double kcl_shomate_coeffs[] = { + 73.59698, 0.0, 0.0, 0.0, 0.0, -443.7341, 175.7209}; + +const double licl_shomate_coeffs[] = { + 73.18025, -9.047232, -0.316390, 0.079587, 0.013594, -417.1314, 157.6711}; }