[Thermo] Make HMWSolution constructible without XML

This commit is contained in:
Ray Speth 2017-08-04 21:27:38 -04:00
parent 4b69c7f265
commit 3790115b99
6 changed files with 598 additions and 1003 deletions

View file

@ -1475,6 +1475,35 @@ public:
* -------------- Utilities -------------------------------
*/
void setBinarySalt(const std::string& sp1, const std::string& sp2,
size_t nParams, double* beta0, double* beta1, double* beta2,
double* Cphi, double alpha1, double alpha2);
void setTheta(const std::string& sp1, const std::string& sp2,
size_t nParams, double* theta);
void setPsi(const std::string& sp1, const std::string& sp2,
const std::string& sp3, size_t nParams, double* psi);
void setLambda(const std::string& sp1, const std::string& sp2,
size_t nParams, double* lambda);
void setMunnn(const std::string& sp, size_t nParams, double* munnn);
void setZeta(const std::string& sp1, const std::string& sp2,
const std::string& sp3, size_t nParams, double* psi);
void setPitzerTempModel(const std::string& model);
void setPitzerRefTemperature(double Tref) {
m_TempPitzerRef = Tref;
}
//! Set the A_Debye parameter. If a negative value is provided, enables
//! calculation of A_Debye using the detailed water equation of state.
void setA_Debye(double A);
void setMaxIonicStrength(double Imax) {
m_maxIionicStrength = Imax;
}
void setCroppingCoefficients(double ln_gamma_k_min, double ln_gamma_k_max,
double ln_gamma_o_min, double ln_gamma_o_max);
virtual void initThermo();
//! Initialize the phase parameters from an XML file.
@ -1759,12 +1788,6 @@ private:
*/
PDSS* m_waterSS;
//! density of standard-state water
/*!
* internal temporary variable
*/
double m_densWaterSS;
//! Pointer to the water property calculator
std::unique_ptr<WaterProps> m_waterProps;
@ -2420,45 +2443,26 @@ private:
*/
void readXMLBinarySalt(XML_Node& BinSalt);
//! Process an XML node called "thetaAnion"
//! Process an XML node called "thetaAnion" or "thetaCation"
/*!
* This node contains all of the parameters necessary to describe the binary
* interactions between two anions.
* interactions between two anions or two cations.
*
* @param BinSalt reference to the XML_Node named thetaAnion containing the
* anion - anion interaction
*/
void readXMLThetaAnion(XML_Node& BinSalt);
void readXMLTheta(XML_Node& BinSalt);
//! Process an XML node called "thetaCation"
/*!
* This node contains all of the parameters necessary to describe the binary
* interactions between two cations.
*
* @param BinSalt reference to the XML_Node named thetaCation containing
* the cation - cation interaction
*/
void readXMLThetaCation(XML_Node& BinSalt);
//! Process an XML node called "psiCommonAnion"
//! Process an XML node called "psiCommonAnion" or "psiCommonCation"
/*!
* This node contains all of the parameters necessary to describe
* the ternary interactions between one anion and two cations.
* the ternary interactions between one anion and two cations or two anions
* and one cation.
*
* @param BinSalt reference to the XML_Node named psiCommonAnion containing
* the anion - cation1 - cation2 interaction
*/
void readXMLPsiCommonAnion(XML_Node& BinSalt);
//! Process an XML node called "psiCommonCation"
/*!
* This node contains all of the parameters necessary to describe
* the ternary interactions between one cation and two anions.
*
* @param BinSalt reference to the XML_Node named psiCommonCation
* containing the cation - anion1 - anion2 interaction
*/
void readXMLPsiCommonCation(XML_Node& BinSalt);
void readXMLPsi(XML_Node& BinSalt);
//! Process an XML node called "lambdaNeutral"
/*!
@ -2491,13 +2495,6 @@ private:
*/
void readXMLZetaCation(const XML_Node& BinSalt);
//! Process an XML node called "croppingCoefficients" for the cropping
//! coefficients values
/*!
* @param acNode Activity Coefficient XML Node
*/
void readXMLCroppingCoefficients(const XML_Node& acNode);
//! Precalculate the IMS Cutoff parameters for typeCutoff = 2
void calcIMSCutoffParams_();

View file

@ -29,10 +29,9 @@ HMWSoln::HMWSoln() :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
m_form_A_Debye(A_DEBYE_WATER),
m_form_A_Debye(A_DEBYE_CONST),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
m_densWaterSS(1000.),
m_molalitiesAreCropped(false),
IMS_X_o_cutoff_(0.2),
IMS_cCut_(0.05),
@ -69,10 +68,9 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
m_form_A_Debye(A_DEBYE_WATER),
m_form_A_Debye(A_DEBYE_CONST),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
m_densWaterSS(1000.),
m_molalitiesAreCropped(false),
IMS_X_o_cutoff_(0.2),
IMS_cCut_(0.05),
@ -106,10 +104,9 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
m_IionicMolality(0.0),
m_maxIionicStrength(100.0),
m_TempPitzerRef(298.15),
m_form_A_Debye(A_DEBYE_WATER),
m_form_A_Debye(A_DEBYE_CONST),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_waterSS(0),
m_densWaterSS(1000.),
m_molalitiesAreCropped(false),
IMS_X_o_cutoff_(0.2),
IMS_cCut_(0.05),
@ -236,10 +233,6 @@ void HMWSoln::calcDensity()
return;
}
// Store the internal density of the water SS. Note, we would have to do
// this for all other species if they had pressure dependent properties.
m_densWaterSS = m_waterSS->density();
// Calculate all of the other standard volumes. Note these are constant for
// now
getPartialMolarVolumes(m_tmpV.data());

File diff suppressed because it is too large Load diff

View file

@ -17,6 +17,7 @@
#include "cantera/thermo/StoichSubstance.h"
#include "cantera/thermo/LatticeSolidPhase.h"
#include "cantera/thermo/IdealSolidSolnPhase.h"
#include "cantera/thermo/HMWSoln.h"
#include "cantera/thermo/NasaPoly2.h"
#include "cantera/thermo/ConstCpPoly.h"
@ -57,11 +58,11 @@ shared_ptr<Species> make_shomate2_species(const std::string& name,
shared_ptr<Species> make_species(const std::string& name,
const std::string& composition, double h298,
double T1, double mu1, double T2, double mu2)
double T1, double mu1, double T2, double mu2, double pref=101325)
{
auto species = make_shared<Species>(name, parseCompString(composition));
double coeffs[] = {2, h298, T1, mu1*GasConstant*T1, T2, mu2*GasConstant*T2};
species->thermo.reset(new Mu0Poly(200, 3500, 101325, coeffs));
species->thermo.reset(new Mu0Poly(200, 3500, pref, coeffs));
return species;
}
@ -550,4 +551,91 @@ TEST(IdealSolidSolnPhase, fromScratch)
EXPECT_NEAR(p.gibbs_mole(), -313642293.1654253, 1e-4);
}
TEST(HMWSoln, fromScratch)
{
// Regression test based on HMW_test_3
HMWSoln p;
p.addUndefinedElements();
auto sH2O = make_species("H2O(l)", "H:2, O:1", h2oliq_nasa_coeffs);
auto sCl = make_species("Cl-", "Cl:1, E:1", 0.0,
298.15, -52.8716, 333.15, -52.8716, 1e5);
sCl->charge = -1;
auto sH = make_species("H+", "H:1, E:-1", 0.0, 298.15, 0.0, 333.15, 0.0, 1e5);
sH->charge = 1;
auto sNa = make_species("Na+", "Na:1, E:-1", 0.0,
298.15, -125.5213, 333.15, -125.5213, 1e5);
sNa->charge = 1;
auto sOH = make_species("OH-", "O:1, H:1, E:1", 0.0,
298.15, -91.523, 333.15, -91.523, 1e5);
sOH->charge = -1;
for (auto& s : {sH2O, sCl, sH, sNa, sOH}) {
p.addSpecies(s);
}
std::unique_ptr<PDSS_Water> ss(new PDSS_Water());
p.installPDSS(0, std::move(ss));
size_t k = 1;
for (double v : {1.3, 1.3, 1.3, 1.3}) {
std::unique_ptr<PDSS_ConstVol> ss(new PDSS_ConstVol());
ss->setMolarVolume(v);
p.installPDSS(k++, std::move(ss));
}
p.setPitzerTempModel("complex");
p.setA_Debye(1.175930);
p.initThermo();
double beta0_nacl[] = {0.0765, 0.008946, -3.3158E-6, -777.03, -4.4706};
double beta1_nacl[] = {0.2664, 6.1608E-5, 1.0715E-6, 0.0, 0.0};
double beta2_nacl[] = {0.0, 0.0, 0.0, 0.0, 0.0};
double cphi_nacl[] = {0.00127, -4.655E-5, 0.0, 33.317, 0.09421};
p.setBinarySalt("Na+", "Cl-", 5, beta0_nacl, beta1_nacl, beta2_nacl,
cphi_nacl, 2.0, 0.0);
double beta0_hcl[] = {0.1775, 0.0, 0.0, 0.0, 0.0};
double beta1_hcl[] = {0.2945, 0.0, 0.0, 0.0, 0.0};
double beta2_hcl[] = {0.0, 0.0, 0.0, 0.0, 0.0};
double cphi_hcl[] = {0.0008, 0.0, 0.0, 0.0, 0.0};
p.setBinarySalt("H+", "Cl-", 5, beta0_hcl, beta1_hcl, beta2_hcl,
cphi_hcl, 2.0, 0.0);
double beta0_naoh[] = {0.0864, 0.0, 0.0, 0.0, 0.0};
double beta1_naoh[] = {0.253, 0.0, 0.0, 0.0, 0.0};
double beta2_naoh[] = {0.0, 0.0, 0.0, 0.0, 0.0};
double cphi_naoh[] = {0.0044, 0.0, 0.0, 0.0, 0.0};
p.setBinarySalt("Na+", "OH-", 5, beta0_naoh, beta1_naoh, beta2_naoh,
cphi_naoh, 2.0, 0.0);
double theta_cloh[] = {-0.05, 0.0, 0.0, 0.0, 0.0};
double psi_nacloh[] = {-0.006, 0.0, 0.0, 0.0, 0.0};
double theta_nah[] = {0.036, 0.0, 0.0, 0.0, 0.0};
double psi_clnah[] = {-0.004, 0.0, 0.0, 0.0, 0.0};
p.setTheta("Cl-", "OH-", 5, theta_cloh);
p.setPsi("Na+", "Cl-", "OH-", 5, psi_nacloh);
p.setTheta("Na+", "H+", 5, theta_nah);
p.setPsi("Cl-", "Na+", "H+", 5, psi_clnah);
p.setMolalitiesByName("Na+:6.0997 Cl-:6.0996986044628 H+:2.1628E-9 OH-:1.3977E-6");
p.setState_TP(150 + 273.15, 101325);
size_t N = p.nSpecies();
vector_fp acMol(N), mf(N), activities(N), moll(N), mu0(N);
p.getMolalityActivityCoefficients(acMol.data());
p.getMoleFractions(mf.data());
p.getActivities(activities.data());
p.getMolalities(moll.data());
p.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);
}
}
} // namespace Cantera

View file

@ -18,8 +18,7 @@
from ThermoPhase that will handle the thermodynamics.
-->
<thermo model="HMW">
<standardConc model="solvent_volume" />
<activityCoefficients model="Pitzer">
<activityCoefficients>
<!-- A_Debye units = sqrt(kg/gmol)
This is adjusted to match the GWB value so
that numerical comparisons can be made
@ -29,8 +28,6 @@
<!-- B_Debye units = sqrt(kg/gmol)/m
-->
<B_Debye> 3.28640E9 </B_Debye>
<ionicRadius default="3.042843" units="Angstroms">
</ionicRadius>
<binarySaltParameters cation="Na+" anion="Cl-">
<beta0> 0.0765 </beta0>
<beta1> 0.2664 </beta1>
@ -60,7 +57,6 @@
</thetaAnion>
<psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
<theta> -0.05 </theta>
<Psi> -0.006 </Psi>
</psiCommonCation>
@ -69,12 +65,10 @@
</thetaCation>
<psiCommonAnion anion="Cl-" cation1="Na+" cation2="H+">
<Theta> 0.036 </Theta>
<Psi> -0.004 </Psi>
</psiCommonAnion>
</activityCoefficients>
<solvent> H2O(L) </solvent>
</thermo>
<elementArray datasrc="elements.xml"> O H C E Fe Si N Na Cl </elementArray>
</phase>

View file

@ -18,8 +18,7 @@
from ThermoPhase that will handle the thermodynamics.
-->
<thermo model="HMW">
<standardConc model="solvent_volume" />
<activityCoefficients model="Pitzer" TempModel="complex1">
<activityCoefficients TempModel="complex1">
<!-- A_Debye units = sqrt(kg/gmol)
This is adjusted to match the GWB value so
that numerical comparisons can be made
@ -28,9 +27,6 @@
<A_Debye> 1.175930 </A_Debye>
<!-- B_Debye units = sqrt(kg/gmol)/m
-->
<B_Debye> 3.28640E9 </B_Debye>
<ionicRadius default="3.042843" units="Angstroms">
</ionicRadius>
<binarySaltParameters cation="Na+" anion="Cl-">
<beta0> 0.0765, 0.008946, -3.3158E-6,
-777.03, -4.4706
@ -44,22 +40,20 @@
</binarySaltParameters>
<binarySaltParameters cation="H+" anion="Cl-">
<beta0> 0.1775, 0.0, 0.0,
0.0, 0.0
<beta0> 0.1775
</beta0>
<beta1> 0.2945, 0.0, 0.0, 0.0, 0.0 </beta1>
<beta2> 0.0, 0.0, 0.0, 0.0, 0.0 </beta2>
<Cphi> 0.0008, 0.0, 0.0,
0.0, 0.0
<beta1> 0.2945 </beta1>
<beta2> 0.0 </beta2>
<Cphi> 0.0008
</Cphi>
<Alpha1> 2.0 </Alpha1>
</binarySaltParameters>
<binarySaltParameters cation="Na+" anion="OH-">
<beta0> 0.0864, 0.0, 0.0, 0.0, 0.0 </beta0>
<beta1> 0.253, 0.0, 0.0, 0.0, 0.0 </beta1>
<beta2> 0.0, 0.0, 0.0, 0.0, 0.0 </beta2>
<Cphi> 0.0044, 0.0, 0.0, 0.0, 0.0 </Cphi>
<beta0> 0.0864 </beta0>
<beta1> 0.253 </beta1>
<beta2> 0.0 </beta2>
<Cphi> 0.0044 </Cphi>
<Alpha1> 2.0 </Alpha1>
</binarySaltParameters>
@ -68,7 +62,6 @@
</thetaAnion>
<psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
<theta> -0.05 </theta>
<Psi> -0.006 </Psi>
</psiCommonCation>
@ -77,12 +70,10 @@
</thetaCation>
<psiCommonAnion anion="Cl-" cation1="Na+" cation2="H+">
<theta> 0.036 </theta>
<Psi> -0.004 </Psi>
</psiCommonAnion>
</activityCoefficients>
<solvent> H2O(L) </solvent>
</thermo>
<elementArray datasrc="elements.xml"> O H C E Fe Si N Na Cl </elementArray>
</phase>