[Thermo] Allow instantiation of MargulesVPSSTP without XML

This commit is contained in:
Ray Speth 2017-06-17 21:57:13 -04:00
parent 4818c87344
commit f69ef44600
4 changed files with 116 additions and 82 deletions

View file

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

View file

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

View file

@ -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<Species> make_species(const std::string& name,
return species;
}
shared_ptr<Species> make_shomate_species(const std::string& name,
const std::string& composition, const double* shomate_coeffs)
{
auto species = make_shared<Species>(name, parseCompString(composition));
species->thermo.reset(new ShomatePoly(200, 3500, 101325, shomate_coeffs));
return species;
}
shared_ptr<Species> 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<PDSS_ConstVol> 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

View file

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