Skip to content

Commit

Permalink
[Thermo] Make PDSS_HKFT constructible without XML
Browse files Browse the repository at this point in the history 
Also fixes PDSS_HKFT to actually make use of the units specified for each input
constant.
  • Loading branch information
@speth
speth committed Aug 14, 2017
1 parent 0ca788b commit 4c630fc
Show file tree
Hide file tree
Showing 5 changed files with 181 additions and 83 deletions.
20 changes: 19 additions & 1 deletion include/cantera/thermo/PDSS_HKFT.h
View file
Expand Up @@ -13,11 +13,11 @@
#define CT_PDSS_HKFT_H

#include "PDSS.h"
#include "WaterProps.h"

namespace Cantera
{
class PDSS_Water;
class WaterProps;

//! Class for pressure dependent standard states corresponding to
//! ionic solutes in electrolyte water.
Expand Down Expand Up @@ -87,6 +87,24 @@ class PDSS_HKFT : public PDSS_Molar
virtual bool useSTITbyPDSS() const { return true; }
virtual void initThermo();

//! Set enthalpy of formation at Pr, Tr [J/kmol]
void setDeltaH0(double dh0);

//! Set Gibbs free energy of formation at Pr, Tr [J/kmol]
void setDeltaG0(double dg0);

//! Set entropy of formation at Pr, Tr [J/kmol/K]
void setS0(double s0);

//! Set "a" coefficients (array of 4 elements). Units of each coefficient
//! are [J/kmol/Pa, J/kmol, J*K/kmol/Pa, J*K/kmol]
void set_a(double* a);

//! Set "c" coefficients (array of 2 elements). Units of each coefficient
//! are [J/kmol/K, J*K/kmol]
void set_c(double* c);
void setOmega(double omega); //!< Set omega [J/kmol]

void setParametersFromXML(const XML_Node& speciesNode);

//! This utility function reports back the type of parameterization and
Expand Down
79 changes: 38 additions & 41 deletions src/thermo/PDSS_HKFT.cpp
View file
Expand Up @@ -316,6 +316,34 @@ void PDSS_HKFT::initThermo()
}
}

void PDSS_HKFT::setDeltaH0(double dh0) {
m_deltaH_formation_tr_pr = dh0 / toSI("cal/gmol");
}

void PDSS_HKFT::setDeltaG0(double dg0) {
m_deltaG_formation_tr_pr = dg0 / toSI("cal/gmol");
}

void PDSS_HKFT::setS0(double s0) {
m_Entrop_tr_pr = s0 / toSI("cal/gmol/K");
}

void PDSS_HKFT::set_a(double* a) {
m_a1 = a[0] / toSI("cal/gmol/bar");
m_a2 = a[1] / toSI("cal/gmol");
m_a3 = a[2] / toSI("cal-K/gmol/bar");
m_a4 = a[3] / toSI("cal-K/gmol");
}

void PDSS_HKFT::set_c(double* c) {
m_c1 = c[0] / toSI("cal/gmol/K");
m_c2 = c[1] / toSI("cal-K/gmol");
}

void PDSS_HKFT::setOmega(double omega) {
m_omega_pr_tr = omega / toSI("cal/gmol");
}

void PDSS_HKFT::setParametersFromXML(const XML_Node& speciesNode)
{
PDSS::setParametersFromXML(speciesNode);
Expand Down Expand Up @@ -357,20 +385,17 @@ void PDSS_HKFT::setParametersFromXML(const XML_Node& speciesNode)
}

if (hh->hasChild("DG0_f_Pr_Tr")) {
doublereal val = getFloat(*hh, "DG0_f_Pr_Tr");
m_deltaG_formation_tr_pr = val;
setDeltaG0(getFloat(*hh, "DG0_f_Pr_Tr", "toSI"));
hasDGO = 1;
}

if (hh->hasChild("DH0_f_Pr_Tr")) {
doublereal val = getFloat(*hh, "DH0_f_Pr_Tr");
m_deltaH_formation_tr_pr = val;
setDeltaH0(getFloat(*hh, "DH0_f_Pr_Tr", "toSI"));
hasDHO = 1;
}

if (hh->hasChild("S0_Pr_Tr")) {
doublereal val = getFloat(*hh, "S0_Pr_Tr");
m_Entrop_tr_pr= val;
setS0(getFloat(*hh, "S0_Pr_Tr", "toSI"));
hasSO = 1;
}

Expand All @@ -384,42 +409,14 @@ void PDSS_HKFT::setParametersFromXML(const XML_Node& speciesNode)
"standardState model for species isn't hkft: "
+ speciesNode.name());
}
if (ss->hasChild("a1")) {
m_a1 = getFloat(*ss, "a1");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a1 field");
}
if (ss->hasChild("a2")) {
m_a2 = getFloat(*ss, "a2");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a2 field");
}
if (ss->hasChild("a3")) {
m_a3 = getFloat(*ss, "a3");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a3 field");
}
if (ss->hasChild("a4")) {
m_a4 = getFloat(*ss, "a4");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a4 field");
}
double a[4] = {getFloat(*ss, "a1", "toSI"), getFloat(*ss, "a2", "toSI"),
getFloat(*ss, "a3", "toSI"), getFloat(*ss, "a4", "toSI")};
set_a(a);

if (ss->hasChild("c1")) {
m_c1 = getFloat(*ss, "c1");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c1 field");
}
if (ss->hasChild("c2")) {
m_c2 = getFloat(*ss, "c2");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c2 field");
}
if (ss->hasChild("omega_Pr_Tr")) {
m_omega_pr_tr = getFloat(*ss, "omega_Pr_Tr");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing omega_Pr_Tr field");
}
double c[2] = {getFloat(*ss, "c1", "toSI"), getFloat(*ss, "c2", "toSI")};
set_c(c);

setOmega(getFloat(*ss, "omega_Pr_Tr", "toSI"));

int isum = hasDGO + hasDHO + hasSO;
if (isum < 2) {
Expand Down
139 changes: 110 additions & 29 deletions test/thermo/phaseConstructors.cpp
View file
Expand Up @@ -18,6 +18,7 @@
#include "cantera/thermo/LatticeSolidPhase.h"
#include "cantera/thermo/IdealSolidSolnPhase.h"
#include "cantera/thermo/HMWSoln.h"
#include "cantera/thermo/PDSS_HKFT.h"

#include "cantera/thermo/NasaPoly2.h"
#include "cantera/thermo/ConstCpPoly.h"
Expand Down Expand Up @@ -551,6 +552,38 @@ TEST(IdealSolidSolnPhase, fromScratch)
EXPECT_NEAR(p.gibbs_mole(), -313642293.1654253, 1e-4);
}

static void set_hmw_interactions(HMWSoln& p) {
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);
}

TEST(HMWSoln, fromScratch)
{
// Regression test based on HMW_test_3
Expand Down Expand Up @@ -583,35 +616,7 @@ TEST(HMWSoln, fromScratch)
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);
set_hmw_interactions(p);
p.setMolalitiesByName("Na+:6.0997 Cl-:6.0996986044628 H+:2.1628E-9 OH-:1.3977E-6");
p.setState_TP(150 + 273.15, 101325);

Expand All @@ -638,4 +643,80 @@ TEST(HMWSoln, fromScratch)
}
}

TEST(HMWSoln, fromScratch_HKFT)
{
HMWSoln p;
p.addUndefinedElements();
auto sH2O = make_species("H2O(l)", "H:2, O:1", h2oliq_nasa_coeffs);
auto sNa = make_species("Na+", "Na:1, E:-1", 0.0,
298.15, -125.5213, 333.15, -125.5213, 1e5);
sNa->charge = 1;
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 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, sNa, sCl, sH, sOH}) {
p.addSpecies(s);
}
double h0[] = {-57433, Undef, 0.0, -54977};
double g0[] = {Undef, -31379, 0.0, -37595};
double s0[] = {13.96, 13.56, Undef, -2.56};
double a[][4] = {{0.1839, -228.5, 3.256, -27260},
{0.4032, 480.1, 5.563, -28470},
{0.0, 0.0, 0.0, 0.0},
{0.12527, 7.38, 1.8423, -27821}};
double c[][2] = {{18.18, -29810}, {-4.4, -57140}, {0.0, 0.0}, {4.15, -103460}};
double omega[] = {33060, 145600, 0.0, 172460};

std::unique_ptr<PDSS_Water> ss(new PDSS_Water());
p.installPDSS(0, std::move(ss));
for (size_t k = 0; k < 4; k++) {
std::unique_ptr<PDSS_HKFT> ss(new PDSS_HKFT());
if (h0[k] != Undef) {
ss->setDeltaH0(h0[k] * toSI("cal/gmol"));
}
if (g0[k] != Undef) {
ss->setDeltaG0(g0[k] * toSI("cal/gmol"));
}
if (s0[k] != Undef) {
ss->setS0(s0[k] * toSI("cal/gmol/K"));
}
a[k][0] *= toSI("cal/gmol/bar");
a[k][1] *= toSI("cal/gmol");
a[k][2] *= toSI("cal-K/gmol/bar");
a[k][3] *= toSI("cal-K/gmol");
c[k][0] *= toSI("cal/gmol/K");
c[k][1] *= toSI("cal-K/gmol");
ss->set_a(a[k]);
ss->set_c(c[k]);
ss->setOmega(omega[k] * toSI("cal/gmol"));
p.installPDSS(k+1, std::move(ss));
}
p.setPitzerTempModel("complex");
p.setA_Debye(-1);
p.initThermo();

set_hmw_interactions(p);
p.setMolalitiesByName("Na+:6.0954 Cl-:6.0954 H+:2.1628E-9 OH-:1.3977E-6");
p.setState_TP(50 + 273.15, 101325);

size_t N = p.nSpecies();
vector_fp mv(N), h(N), mu(N), ac(N), acoeff(N);
p.getPartialMolarVolumes(mv.data());
p.getPartialMolarEnthalpies(h.data());
p.getChemPotentials(mu.data());
p.getActivities(ac.data());
p.getActivityCoefficients(acoeff.data());

double mvRef[] = {0.01815224, 0.00157182, 0.01954605, 0.00173137, -0.0020266};

for (size_t k = 0; k < N; k++) {
EXPECT_NEAR(mv[k], mvRef[k], 2e-8);
}
}

} // namespace Cantera
2 changes: 2 additions & 0 deletions test/thermo/thermoParameterizations.cpp
View file
Expand Up @@ -4,8 +4,10 @@
#include "cantera/thermo/ConstCpPoly.h"
#include "cantera/thermo/NasaPoly2.h"
#include "cantera/thermo/ShomatePoly.h"
#include "cantera/thermo/PDSS_HKFT.h"
#include "cantera/base/stringUtils.h"
#include "thermo_data.h"
#include <sstream>

using namespace Cantera;

Expand Down
24 changes: 12 additions & 12 deletions test_problems/simpleTransport/HMW_NaCl_pdss.xml
View file
Expand Up @@ -141,10 +141,10 @@
<standardState model="HKFT">
<a1 units="cal/gmol/bar"> 0.1839 </a1>
<a2 units="cal/gmol"> -228.5 </a2>
<a3 units="cal K/gmol/bar"> 3.256 </a3>
<a4 units="cal K/gmol"> -27260. </a4>
<a3 units="cal-K/gmol/bar"> 3.256 </a3>
<a4 units="cal-K/gmol"> -27260. </a4>
<c1 units="cal/gmol/K"> 18.18 </c1>
<c2 units="cal K/gmol"> -29810. </c2>
<c2 units="cal-K/gmol"> -29810. </c2>
<omega_Pr_Tr units="cal/gmol"> 33060. </omega_Pr_Tr>
</standardState>
<transport>
Expand All @@ -168,10 +168,10 @@
<standardState model="HKFT">
<a1 units="cal/gmol/bar"> 0.4032 </a1>
<a2 units="cal/gmol"> 480.1 </a2>
<a3 units="cal K/gmol/bar"> 5.563 </a3>
<a4 units="cal K/gmol"> -28470. </a4>
<a3 units="cal-K/gmol/bar"> 5.563 </a3>
<a4 units="cal-K/gmol"> -28470. </a4>
<c1 units="cal/gmol/K"> -4.4 </c1>
<c2 units="cal K/gmol"> -57140. </c2>
<c2 units="cal-K/gmol"> -57140. </c2>
<omega_Pr_Tr units="cal/gmol"> 145600. </omega_Pr_Tr>
</standardState>
<transport>
Expand All @@ -195,10 +195,10 @@
<standardState model="HKFT">
<a1 units="cal/gmol/bar"> 0.0 </a1>
<a2 units="cal/gmol"> 0.0 </a2>
<a3 units="cal K/gmol/bar"> 0.0 </a3>
<a4 units="cal K/gmol"> 0.0 </a4>
<a3 units="cal-K/gmol/bar"> 0.0 </a3>
<a4 units="cal-K/gmol"> 0.0 </a4>
<c1 units="cal/gmol/K"> 0.0 </c1>
<c2 units="cal K/gmol"> 0.0 </c2>
<c2 units="cal-K/gmol"> 0.0 </c2>
<omega_Pr_Tr units="cal/gmol"> 0.0 </omega_Pr_Tr>
</standardState>
<transport>
Expand All @@ -223,10 +223,10 @@
<standardState model="HKFT">
<a1 units="cal/gmol/bar"> 0.12527 </a1>
<a2 units="cal/gmol"> 7.38 </a2>
<a3 units="cal K/gmol/bar"> 1.8423 </a3>
<a4 units="cal K/gmol"> -27821 </a4>
<a3 units="cal-K/gmol/bar"> 1.8423 </a3>
<a4 units="cal-K/gmol"> -27821 </a4>
<c1 units="cal/gmol/K"> 4.15 </c1>
<c2 units="cal K/gmol"> -103460. </c2>
<c2 units="cal-K/gmol"> -103460. </c2>
<omega_Pr_Tr units="cal/gmol"> 172460. </omega_Pr_Tr>
</standardState>
<transport>
Expand Down

0 comments on commit 4c630fc

Please sign in to comment.