31 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
32 m_IionicMolality(0.0),
33 m_maxIionicStrength(100.0),
34 m_TempPitzerRef(298.15),
35 m_form_A_Debye(A_DEBYE_CONST),
38 m_molalitiesAreCropped(false),
56 CROP_ln_gamma_o_min(-6.0),
57 CROP_ln_gamma_o_max(3.0),
58 CROP_ln_gamma_k_min(-5.0),
59 CROP_ln_gamma_k_max(15.0),
69 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
70 m_IionicMolality(0.0),
71 m_maxIionicStrength(100.0),
72 m_TempPitzerRef(298.15),
73 m_form_A_Debye(A_DEBYE_CONST),
76 m_molalitiesAreCropped(false),
94 CROP_ln_gamma_o_min(-6.0),
95 CROP_ln_gamma_o_max(3.0),
96 CROP_ln_gamma_k_min(-5.0),
97 CROP_ln_gamma_k_max(15.0),
104 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
105 m_IionicMolality(0.0),
106 m_maxIionicStrength(100.0),
107 m_TempPitzerRef(298.15),
108 m_form_A_Debye(A_DEBYE_CONST),
111 m_molalitiesAreCropped(false),
112 IMS_X_o_cutoff_(0.2),
129 CROP_ln_gamma_o_min(-6.0),
130 CROP_ln_gamma_o_max(3.0),
131 CROP_ln_gamma_k_min(-5.0),
132 CROP_ln_gamma_k_max(15.0),
151 for (
size_t k = 0; k <
m_kk; k++) {
163 size_t kcation =
npos;
164 double xcation = 0.0;
165 size_t kanion =
npos;
166 for (
size_t k = 0; k <
m_kk; k++) {
172 }
else if (
charge(k) < 0.0) {
173 if (
m_tmpV[k] > xcation) {
179 if (kcation ==
npos || kanion ==
npos) {
182 double xuse = xcation;
184 if (xanion < xcation) {
186 if (
charge(kcation) != 1.0) {
190 if (
charge(kanion) != 1.0) {
194 xuse = xuse / factor;
223 return cp - beta * beta * tt * molarV / kappa_t;
246 if (rho != dens_old) {
248 "Density is not an independent variable");
255 "Density is not an independent variable");
267 for (
size_t k = 1; k <
m_kk; k++) {
276 double mvSolvent =
m_tmpV[0];
280 return 1.0 / mvSolvent;
289 s_update_lnMolalityActCoeff();
292 for (
size_t k = 1; k <
m_kk; k++) {
303 s_update_lnMolalityActCoeff();
305 for (
size_t k = 0; k <
m_kk; k++) {
306 acMolality[k] = exp(acMolality[k]);
322 s_update_lnMolalityActCoeff();
324 for (
size_t k = 1; k <
m_kk; k++) {
338 for (
size_t k = 0; k <
m_kk; k++) {
344 s_update_lnMolalityActCoeff();
346 for (
size_t k = 0; k <
m_kk; k++) {
358 for (
size_t k = 0; k <
m_kk; k++) {
364 s_update_lnMolalityActCoeff();
369 for (
size_t k = 1; k <
m_kk; k++) {
381 for (
size_t k = 0; k <
m_kk; k++) {
392 s_update_lnMolalityActCoeff();
394 for (
size_t k = 0; k <
m_kk; k++) {
402 for (
size_t k = 0; k <
m_kk; k++) {
408 s_update_lnMolalityActCoeff();
411 for (
size_t k = 0; k <
m_kk; k++) {
429 static void check_nParams(
const std::string& method,
size_t nParams,
430 size_t m_formPitzerTemp)
432 if (m_formPitzerTemp == PITZER_TEMP_CONSTANT && nParams != 1) {
433 throw CanteraError(method,
"'constant' temperature model requires one" 434 " coefficient for each of parameter, but {} were given", nParams);
435 }
else if (m_formPitzerTemp == PITZER_TEMP_LINEAR && nParams != 2) {
436 throw CanteraError(method,
"'linear' temperature model requires two" 437 " coefficients for each parameter, but {} were given", nParams);
439 if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1 && nParams != 5) {
440 throw CanteraError(method,
"'complex' temperature model requires five" 441 " coefficients for each parameter, but {} were given", nParams);
445 void HMWSoln::setBinarySalt(
const std::string& sp1,
const std::string& sp2,
446 size_t nParams,
double* beta0,
double* beta1,
double* beta2,
447 double* Cphi,
double alpha1,
double alpha2)
452 throw CanteraError(
"HMWSoln::setBinarySalt",
"Species '{}' not found", sp1);
453 }
else if (k2 ==
npos) {
454 throw CanteraError(
"HMWSoln::setBinarySalt",
"Species '{}' not found", sp2);
459 throw CanteraError(
"HMWSoln::setBinarySalt",
"Species '{}' and '{}' " 460 "do not have opposite charges ({}, {})", sp1, sp2,
470 for (
size_t n = 0; n < nParams; n++) {
476 m_Alpha1MX_ij[c] = alpha1;
480 void HMWSoln::setTheta(
const std::string& sp1,
const std::string& sp2,
481 size_t nParams,
double* theta)
486 throw CanteraError(
"HMWSoln::setTheta",
"Species '{}' not found", sp1);
487 }
else if (k2 ==
npos) {
488 throw CanteraError(
"HMWSoln::setTheta",
"Species '{}' not found", sp2);
491 throw CanteraError(
"HMWSoln::setTheta",
"Species '{}' and '{}' " 492 "should both have the same (non-zero) charge ({}, {})", sp1, sp2,
498 for (
size_t n = 0; n < nParams; n++) {
503 void HMWSoln::setPsi(
const std::string& sp1,
const std::string& sp2,
504 const std::string& sp3,
size_t nParams,
double* psi)
510 throw CanteraError(
"HMWSoln::setPsi",
"Species '{}' not found", sp1);
511 }
else if (k2 ==
npos) {
512 throw CanteraError(
"HMWSoln::setPsi",
"Species '{}' not found", sp2);
513 }
else if (k3 ==
npos) {
514 throw CanteraError(
"HMWSoln::setPsi",
"Species '{}' not found", sp3);
519 throw CanteraError(
"HMWSoln::setPsi",
"All species must be ions and" 520 " must include at least one cation and one anion, but given species" 521 " (charges) were: {} ({}), {} ({}), and {} ({}).",
532 for (
size_t n = 0; n < nParams; n++) {
539 void HMWSoln::setLambda(
const std::string& sp1,
const std::string& sp2,
540 size_t nParams,
double* lambda)
545 throw CanteraError(
"HMWSoln::setLambda",
"Species '{}' not found", sp1);
546 }
else if (k2 ==
npos) {
547 throw CanteraError(
"HMWSoln::setLambda",
"Species '{}' not found", sp2);
551 throw CanteraError(
"HMWSoln::setLambda",
"Expected at least one neutral" 552 " species, but given species (charges) were: {} ({}) and {} ({}).",
559 size_t c = k1*
m_kk + k2;
560 for (
size_t n = 0; n < nParams; n++) {
566 void HMWSoln::setMunnn(
const std::string& sp,
size_t nParams,
double* munnn)
570 throw CanteraError(
"HMWSoln::setMunnn",
"Species '{}' not found", sp);
574 throw CanteraError(
"HMWSoln::setMunnn",
"Expected a neutral species," 575 " got {} ({}).", sp,
charge(k));
578 for (
size_t n = 0; n < nParams; n++) {
584 void HMWSoln::setZeta(
const std::string& sp1,
const std::string& sp2,
585 const std::string& sp3,
size_t nParams,
double* psi)
591 throw CanteraError(
"HMWSoln::setZeta",
"Species '{}' not found", sp1);
592 }
else if (k2 ==
npos) {
593 throw CanteraError(
"HMWSoln::setZeta",
"Species '{}' not found", sp2);
594 }
else if (k3 ==
npos) {
595 throw CanteraError(
"HMWSoln::setZeta",
"Species '{}' not found", sp3);
600 throw CanteraError(
"HMWSoln::setZeta",
"Requires one neutral species, " 601 "one cation, and one anion, but given species (charges) were: " 602 "{} ({}), {} ({}), and {} ({}).",
609 }
else if (
charge(k3) == 0) {
621 for (
size_t n = 0; n < nParams; n++) {
627 void HMWSoln::setPitzerTempModel(
const std::string& model)
636 throw CanteraError(
"HMWSoln::setPitzerTempModel",
637 "Unknown Pitzer ActivityCoeff Temp model: {}", model);
651 void HMWSoln::setCroppingCoefficients(
double ln_gamma_k_min,
652 double ln_gamma_k_max,
double ln_gamma_o_min,
double ln_gamma_o_max)
654 CROP_ln_gamma_k_min = ln_gamma_k_min;
655 CROP_ln_gamma_k_max = ln_gamma_k_max;
656 CROP_ln_gamma_o_min = ln_gamma_o_min;
657 CROP_ln_gamma_o_max = ln_gamma_o_max;
665 for (
int i = 0; i < 17; i++) {
687 for (
size_t k = 0; k <
m_kk; k++) {
689 if (fabs(mf[k] *
charge(k)) > MaxC) {
696 if (fabs(sum) > 1.0E-30) {
698 if (mf[kHp] > sum * 1.1) {
712 if (mf[kOHm] > -sum * 1.1) {
724 if (notDone && kMaxC !=
npos) {
725 if (mf[kMaxC] > (1.1 * sum /
charge(kMaxC))) {
726 mf[kMaxC] -= sum /
charge(kMaxC);
727 mf[0] += sum /
charge(kMaxC);
748 if (id_.size() > 0) {
749 string idp = phaseNode.
id();
752 "phasenode and Id are incompatible");
757 if (!phaseNode.
hasChild(
"thermo")) {
759 "no thermo XML node");
765 if (thermoNode.
hasChild(
"activityCoefficients")) {
770 string formString = scNode.
attrib(
"TempModel");
771 if (formString !=
"") {
772 setPitzerTempModel(formString);
778 formString = scNode.
attrib(
"TempReference");
779 if (formString !=
"") {
790 if (thermoNode.
hasChild(
"activityCoefficients")) {
804 if (acNode.
hasChild(
"maxIonicStrength")) {
805 setMaxIonicStrength(
getFloat(acNode,
"maxIonicStrength"));
808 for (
const auto& xmlACChild : acNode.
children()) {
809 string nodeName = xmlACChild->name();
832 if (acNode.
hasChild(
"croppingCoefficients")) {
834 setCroppingCoefficients(
835 getFloat(cropNode.
child(
"ln_gamma_k_min"),
"pureSolventValue"),
836 getFloat(cropNode.
child(
"ln_gamma_k_max"),
"pureSolventValue"),
837 getFloat(cropNode.
child(
"ln_gamma_o_min"),
"pureSolventValue"),
838 getFloat(cropNode.
child(
"ln_gamma_o_max"),
"pureSolventValue"));
849 if (tempArg != -1.0) {
853 if (presArg != -1.0) {
872 throw CanteraError(
"HMWSoln::A_Debye_TP",
"shouldn't be here");
880 if (tempArg != -1.0) {
884 if (presArg != -1.0) {
896 throw CanteraError(
"HMWSoln::dA_DebyedT_TP",
"shouldn't be here");
904 if (tempArg != -1.0) {
908 if (presArg != -1.0) {
928 throw CanteraError(
"HMWSoln::dA_DebyedP_TP",
"shouldn't be here");
936 double dAphidT = dAdT /3.0;
938 if (tempArg != -1.0) {
947 double dAphidP = dAdP /3.0;
949 if (tempArg != -1.0) {
958 if (tempArg != -1.0) {
963 double d2Aphi = d2 / 3.0;
964 return 2.0 * A_L / T + 4.0 *
GasConstant * T * T *d2Aphi;
970 if (tempArg != -1.0) {
974 if (presArg != -1.0) {
986 throw CanteraError(
"HMWSoln::d2A_DebyedT2_TP",
"shouldn't be here");
1000 size_t maxCounterIJlen = 1 + (
m_kk-1) * (
m_kk-2) / 2;
1003 int TCoeffLength = 1;
1034 m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
1075 m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1087 m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1096 m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1108 void HMWSoln::s_update_lnMolalityActCoeff()
const 1136 double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1137 double lnxs = log(xx);
1139 for (
size_t k = 1; k <
m_kk; k++) {
1174 doublereal Imax = 0.0;
1177 for (
size_t k = 0; k <
m_kk; k++) {
1183 if (cropMethod == 0) {
1190 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1191 double charge_i =
charge(i);
1192 double abs_charge_i = fabs(charge_i);
1193 if (charge_i == 0.0) {
1196 for (
size_t j = (i+1); j <
m_kk; j++) {
1197 double charge_j =
charge(j);
1198 double abs_charge_j = fabs(charge_j);
1201 if (charge_i * charge_j < 0) {
1206 if (Imax > Iac_max) {
1210 if (Imax > Iac_max) {
1215 if (Imax > Iac_max) {
1219 if (Imax > Iac_max) {
1229 for (
int times = 0; times< 10; times++) {
1230 double anion_charge = 0.0;
1231 double cation_charge = 0.0;
1232 size_t anion_contrib_max_i =
npos;
1233 double anion_contrib_max = -1.0;
1234 size_t cation_contrib_max_i =
npos;
1235 double cation_contrib_max = -1.0;
1236 for (
size_t i = 0; i <
m_kk; i++) {
1237 double charge_i =
charge(i);
1238 if (charge_i < 0.0) {
1240 anion_charge += anion_contrib;
1241 if (anion_contrib > anion_contrib_max) {
1242 anion_contrib_max = anion_contrib;
1243 anion_contrib_max_i = i;
1245 }
else if (charge_i > 0.0) {
1247 cation_charge += cation_contrib;
1248 if (cation_contrib > cation_contrib_max) {
1249 cation_contrib_max = cation_contrib;
1250 cation_contrib_max_i = i;
1254 double total_charge = cation_charge - anion_charge;
1255 if (total_charge > 1.0E-8) {
1256 double desiredCrop = total_charge/
charge(cation_contrib_max_i);
1258 if (desiredCrop < maxCrop) {
1264 }
else if (total_charge < -1.0E-8) {
1265 double desiredCrop = total_charge/
charge(anion_contrib_max_i);
1267 if (desiredCrop < maxCrop) {
1279 if (cropMethod == 1) {
1282 double xmolSolvent = molF[0];
1288 double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1289 double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1291 for (
size_t k = 0; k <
m_kk; k++) {
1299 for (
size_t k = 0; k <
m_kk; k++) {
1304 for (
size_t k = 0; k <
m_kk; k++) {
1317 for (
size_t i = 0; i <
m_kk; i++) {
1319 size_t nc =
m_kk * i;
1323 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1324 size_t n =
m_kk * i + i;
1326 for (
size_t j = (i+1); j <
m_kk; j++) {
1328 size_t nc =
m_kk * i + j;
1338 if (BinSalt.
name() !=
"binarySaltParameters") {
1340 "Incorrect name for processing this routine: " + BinSalt.
name());
1343 string iName = BinSalt.
attrib(
"cation");
1345 throw CanteraError(
"HMWSoln::readXMLBinarySalt",
"no cation attrib");
1347 string jName = BinSalt.
attrib(
"anion");
1349 throw CanteraError(
"HMWSoln::readXMLBinarySalt",
"no anion attrib");
1363 if (beta0.size() != beta1.size() || beta0.size() != beta2.size() ||
1364 beta0.size() != Cphi.size()) {
1365 throw CanteraError(
"HMWSoln::readXMLBinarySalt",
"Inconsistent" 1366 " array sizes ({}, {}, {}, {})", beta0.size(), beta1.size(),
1367 beta2.size(), Cphi.size());
1370 beta0.resize(5, 0.0);
1371 beta1.resize(5, 0.0);
1372 beta2.resize(5, 0.0);
1373 Cphi.resize(5, 0.0);
1375 double alpha1 =
getFloat(BinSalt,
"Alpha1");
1376 double alpha2 = 0.0;
1378 setBinarySalt(iName, jName, beta0.size(), beta0.data(), beta1.data(),
1379 beta2.data(), Cphi.data(), alpha1, alpha2);
1384 string ispName, jspName;
1385 if (node.
name() ==
"thetaAnion") {
1386 ispName = node.
attrib(
"anion1");
1387 if (ispName ==
"") {
1388 throw CanteraError(
"HMWSoln::readXMLTheta",
"no anion1 attrib");
1390 jspName = node.
attrib(
"anion2");
1391 if (jspName ==
"") {
1392 throw CanteraError(
"HMWSoln::readXMLTheta",
"no anion2 attrib");
1394 }
else if (node.
name() ==
"thetaCation") {
1395 ispName = node.
attrib(
"cation1");
1396 if (ispName ==
"") {
1397 throw CanteraError(
"HMWSoln::readXMLTheta",
"no cation1 attrib");
1399 jspName = node.
attrib(
"cation2");
1400 if (jspName ==
"") {
1401 throw CanteraError(
"HMWSoln::readXMLTheta",
"no cation2 attrib");
1405 "Incorrect name for processing this routine: " + node.
name());
1417 theta.resize(5, 0.0);
1419 setTheta(ispName, jspName, theta.size(), theta.data());
1424 string iName, jName, kName;
1425 if (node.
name() ==
"psiCommonCation") {
1426 kName = node.
attrib(
"cation");
1428 throw CanteraError(
"HMWSoln::readXMLPsi",
"no cation attrib");
1430 iName = node.
attrib(
"anion1");
1432 throw CanteraError(
"HMWSoln::readXMLPsi",
"no anion1 attrib");
1434 jName = node.
attrib(
"anion2");
1436 throw CanteraError(
"HMWSoln::readXMLPsi",
"no anion2 attrib");
1438 }
else if (node.
name() ==
"psiCommonAnion") {
1439 kName = node.
attrib(
"anion");
1441 throw CanteraError(
"HMWSoln::readXMLPsi",
"no anion attrib");
1443 iName = node.
attrib(
"cation1");
1445 throw CanteraError(
"HMWSoln::readXMLPsi",
"no cation1 attrib");
1447 jName = node.
attrib(
"cation2");
1449 throw CanteraError(
"HMWSoln::readXMLPsi",
"no cation2 attrib");
1453 "Incorrect name for processing this routine: " + node.
name());
1468 setPsi(iName, jName, kName, psi.size(), psi.data());
1474 if (node.
name() !=
"lambdaNeutral") {
1476 "Incorrect name for processing this routine: " + node.
name());
1478 string iName = node.
attrib(
"species1");
1480 throw CanteraError(
"HMWSoln::readXMLLambdaNeutral",
"no species1 attrib");
1482 string jName = node.
attrib(
"species2");
1484 throw CanteraError(
"HMWSoln::readXMLLambdaNeutral",
"no species2 attrib");
1496 lambda.resize(5, 0.0);
1498 setLambda(iName, jName, lambda.size(), lambda.data());
1503 if (node.
name() !=
"MunnnNeutral") {
1505 "Incorrect name for processing this routine: " + node.
name());
1507 string iName = node.
attrib(
"species1");
1509 throw CanteraError(
"HMWSoln::readXMLMunnnNeutral",
"no species1 attrib");
1521 munnn.resize(5, 0.0);
1523 setMunnn(iName, munnn.size(), munnn.data());
1528 if (node.
name() !=
"zetaCation") {
1530 "Incorrect name for processing this routine: " + node.
name());
1533 string iName = node.
attrib(
"neutral");
1535 throw CanteraError(
"HMWSoln::readXMLZetaCation",
"no neutral attrib");
1537 string jName = node.
attrib(
"cation1");
1539 throw CanteraError(
"HMWSoln::readXMLZetaCation",
"no cation1 attrib");
1541 string kName = node.
attrib(
"anion1");
1543 throw CanteraError(
"HMWSoln::readXMLZetaCation",
"no anion1 attrib");
1556 zeta.resize(5, 0.0);
1558 setZeta(iName, jName, kName, zeta.size(), zeta.data());
1563 double IMS_gamma_o_min_ = 1.0E-5;
1564 double IMS_gamma_k_min_ = 10.0;
1565 double IMS_slopefCut_ = 0.6;
1567 IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
1569 bool converged =
false;
1571 for (
int its = 0; its < 100 && !converged; its++) {
1573 IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) -IMS_efCut_;
1574 IMS_bfCut_ = IMS_afCut_ /
IMS_cCut_ + IMS_slopefCut_ - 1.0;
1580 IMS_efCut_ = - eterm * tmp;
1581 if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
1587 " failed to converge on the f polynomial");
1590 double f_0 = IMS_afCut_ + IMS_efCut_;
1591 double f_prime_0 = 1.0 - IMS_afCut_ /
IMS_cCut_ + IMS_bfCut_;
1593 for (
int its = 0; its < 100 && !converged; its++) {
1595 double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
1596 IMS_agCut_ = exp(lng_0) - IMS_egCut_;
1603 IMS_egCut_ = - eterm * tmp;
1604 if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
1610 " failed to converge on the g polynomial");
1616 double MC_X_o_min_ = 0.35;
1618 double MC_slopepCut_ = 0.02;
1622 MC_apCut_ = MC_X_o_min_;
1624 bool converged =
false;
1627 for (
int its = 0; its < 500 && !converged; its++) {
1629 MC_apCut_ = damp *(MC_X_o_min_ - MC_epCut_) + (1-damp) * MC_apCut_;
1630 double MC_bpCutNew = MC_apCut_ / MC_cpCut_ + MC_slopepCut_ - 1.0;
1631 MC_bpCut_ = damp * MC_bpCutNew + (1-damp) * MC_bpCut_;
1632 double MC_dpCutNew = ((- MC_apCut_/MC_cpCut_ + MC_bpCut_ - MC_bpCut_ *
MC_X_o_cutoff_/MC_cpCut_)
1635 MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
1638 MC_epCut_ = - eterm * tmp;
1639 double diff = MC_epCut_ - oldV;
1640 if (fabs(diff) < 1.0E-14) {
1646 " failed to converge on the p polynomial");
1653 const double twoT = 2.0 * T;
1654 const double invT = 1.0 / T;
1655 const double invT2 = invT * invT;
1656 const double twoinvT3 = 2.0 * invT * invT2;
1657 double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
1667 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1668 for (
size_t j = (i+1); j <
m_kk; j++) {
1670 size_t n =
m_kk*i + j;
1680 case PITZER_TEMP_CONSTANT:
1682 case PITZER_TEMP_LINEAR:
1685 + beta0MX_coeff[1]*tlin;
1689 + beta1MX_coeff[1]*tlin;
1693 + beta2MX_coeff[1]*tlin;
1697 + CphiMX_coeff[1]*tlin;
1700 m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1705 case PITZER_TEMP_COMPLEX1:
1707 + beta0MX_coeff[1]*tlin
1708 + beta0MX_coeff[2]*tquad
1709 + beta0MX_coeff[3]*tinv
1710 + beta0MX_coeff[4]*tln;
1712 + beta1MX_coeff[1]*tlin
1713 + beta1MX_coeff[2]*tquad
1714 + beta1MX_coeff[3]*tinv
1715 + beta1MX_coeff[4]*tln;
1717 + beta2MX_coeff[1]*tlin
1718 + beta2MX_coeff[2]*tquad
1719 + beta2MX_coeff[3]*tinv
1720 + beta2MX_coeff[4]*tln;
1722 + CphiMX_coeff[1]*tlin
1723 + CphiMX_coeff[2]*tquad
1724 + CphiMX_coeff[3]*tinv
1725 + CphiMX_coeff[4]*tln;
1727 + Theta_coeff[1]*tlin
1728 + Theta_coeff[2]*tquad
1729 + Theta_coeff[3]*tinv
1730 + Theta_coeff[4]*tln;
1732 + beta0MX_coeff[2]*twoT
1733 - beta0MX_coeff[3]*invT2
1734 + beta0MX_coeff[4]*invT;
1736 + beta1MX_coeff[2]*twoT
1737 - beta1MX_coeff[3]*invT2
1738 + beta1MX_coeff[4]*invT;
1740 + beta2MX_coeff[2]*twoT
1741 - beta2MX_coeff[3]*invT2
1742 + beta2MX_coeff[4]*invT;
1744 + CphiMX_coeff[2]*twoT
1745 - CphiMX_coeff[3]*invT2
1746 + CphiMX_coeff[4]*invT;
1748 + Theta_coeff[2]*twoT
1749 - Theta_coeff[3]*invT2
1750 + Theta_coeff[4]*invT;
1754 + beta0MX_coeff[2]*2.0
1755 + beta0MX_coeff[3]*twoinvT3
1756 - beta0MX_coeff[4]*invT2;
1758 + beta1MX_coeff[2]*2.0
1759 + beta1MX_coeff[3]*twoinvT3
1760 - beta1MX_coeff[4]*invT2;
1762 + beta2MX_coeff[2]*2.0
1763 + beta2MX_coeff[3]*twoinvT3
1764 - beta2MX_coeff[4]*invT2;
1766 + CphiMX_coeff[2]*2.0
1767 + CphiMX_coeff[3]*twoinvT3
1768 - CphiMX_coeff[4]*invT2;
1770 + Theta_coeff[2]*2.0
1771 + Theta_coeff[3]*twoinvT3
1772 - Theta_coeff[4]*invT2;
1782 for (
size_t i = 1; i <
m_kk; i++) {
1784 for (
size_t j = 1; j <
m_kk; j++) {
1785 size_t n = i *
m_kk + j;
1788 case PITZER_TEMP_CONSTANT:
1791 case PITZER_TEMP_LINEAR:
1792 m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1796 case PITZER_TEMP_COMPLEX1:
1798 + Lambda_coeff[1]*tlin
1799 + Lambda_coeff[2]*tquad
1800 + Lambda_coeff[3]*tinv
1801 + Lambda_coeff[4]*tln;
1804 + Lambda_coeff[2]*twoT
1805 - Lambda_coeff[3]*invT2
1806 + Lambda_coeff[4]*invT;
1810 + Lambda_coeff[3]*twoinvT3
1811 - Lambda_coeff[4]*invT2;
1817 case PITZER_TEMP_CONSTANT:
1820 case PITZER_TEMP_LINEAR:
1821 m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1825 case PITZER_TEMP_COMPLEX1:
1837 + Mu_coeff[3]*twoinvT3
1838 - Mu_coeff[4]*invT2;
1846 case PITZER_TEMP_CONSTANT:
1847 for (
size_t i = 1; i <
m_kk; i++) {
1848 for (
size_t j = 1; j <
m_kk; j++) {
1849 for (
size_t k = 1; k <
m_kk; k++) {
1857 case PITZER_TEMP_LINEAR:
1858 for (
size_t i = 1; i <
m_kk; i++) {
1859 for (
size_t j = 1; j <
m_kk; j++) {
1860 for (
size_t k = 1; k <
m_kk; k++) {
1863 m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1870 case PITZER_TEMP_COMPLEX1:
1871 for (
size_t i = 1; i <
m_kk; i++) {
1872 for (
size_t j = 1; j <
m_kk; j++) {
1873 for (
size_t k = 1; k <
m_kk; k++) {
1878 + Psi_coeff[2]*tquad
1883 - Psi_coeff[3]*invT2
1884 + Psi_coeff[4]*invT;
1887 + Psi_coeff[3]*twoinvT3
1888 - Psi_coeff[4]*invT2;
1906 double etheta[5][5], etheta_prime[5][5], sqrtIs;
1913 double molarcharge = 0.0;
1917 double molalitysumUncropped = 0.0;
1923 for (
size_t n = 1; n <
m_kk; n++) {
1927 molarcharge += fabs(
charge(n)) * molality[n];
1942 for (
int z1 = 1; z1 <=4; z1++) {
1943 for (
int z2 =1; z2 <=4; z2++) {
1944 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
1951 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1952 for (
size_t j = (i+1); j <
m_kk; j++) {
1954 size_t n =
m_kk*i + j;
1960 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
1961 if (x1 > 1.0E-100) {
1962 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
1964 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
1972 if (x2 > 1.0E-100) {
1973 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
1975 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
1990 for (
size_t i = 1; i <
m_kk - 1; i++) {
1991 for (
size_t j = i+1; j <
m_kk; j++) {
1993 size_t n =
m_kk*i + j;
2003 if (Is > 1.0E-150) {
2020 for (
size_t i = 1; i <
m_kk-1; i++) {
2021 for (
size_t j = i+1; j <
m_kk; j++) {
2023 size_t n =
m_kk*i + j;
2039 for (
size_t i = 1; i <
m_kk-1; i++) {
2040 for (
size_t j = i+1; j <
m_kk; j++) {
2042 size_t n =
m_kk*i + j;
2048 int z1 = (int) fabs(
charge(i));
2049 int z2 = (int) fabs(
charge(j));
2064 double F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2065 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2066 for (
size_t i = 1; i <
m_kk-1; i++) {
2067 for (
size_t j = i+1; j <
m_kk; j++) {
2069 size_t n =
m_kk*i + j;
2086 for (
size_t i = 1; i <
m_kk; i++) {
2099 for (
size_t j = 1; j <
m_kk; j++) {
2101 size_t n =
m_kk*i + j;
2106 sum1 += molality[j] *
2112 for (
size_t k = j+1; k <
m_kk; k++) {
2116 sum3 += molality[j]*molality[k]*
m_Psi_ijk[n];
2125 sum2 += molality[j]*(2.0*
m_Phi_IJ[counterIJ]);
2127 for (
size_t k = 1; k <
m_kk; k++) {
2131 sum2 += molality[j]*molality[k]*
m_Psi_ijk[n];
2136 sum4 += (fabs(
charge(i))*
2137 molality[j]*molality[k]*
m_CMX_IJ[counterIJ2]);
2147 for (
size_t k = 1; k <
m_kk; k++) {
2154 sum5 += molality[j]*molality[k]*zeta;
2178 for (
size_t j = 1; j <
m_kk; j++) {
2180 size_t n =
m_kk*i + j;
2185 sum1 += molality[j]*
2188 for (
size_t k = j+1; k <
m_kk; k++) {
2192 sum3 += molality[j]*molality[k]*
m_Psi_ijk[n];
2202 sum2 += molality[j]*(2.0*
m_Phi_IJ[counterIJ]);
2204 for (
size_t k = 1; k <
m_kk; k++) {
2208 sum2 += molality[j]*molality[k]*
m_Psi_ijk[n];
2213 molality[j]*molality[k]*
m_CMX_IJ[counterIJ2];
2222 for (
size_t k = 1; k <
m_kk; k++) {
2230 sum5 += molality[j]*molality[k]*zeta;
2246 for (
size_t j = 1; j <
m_kk; j++) {
2250 for (
size_t k = 1; k <
m_kk; k++) {
2253 sum3 += molality[j]*molality[k]*
m_Psi_ijk[n];
2258 double sum2 = 3.0 * molality[i]* molality[i] *
m_Mu_nnn[i];
2280 double term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2282 for (
size_t j = 1; j <
m_kk; j++) {
2285 for (
size_t k = 1; k <
m_kk; k++) {
2288 size_t n =
m_kk*j + k;
2291 sum1 += molality[j]*molality[k]*
2296 for (
size_t k = j+1; k <
m_kk; k++) {
2297 if (j == (
m_kk-1)) {
2299 throw CanteraError(
"HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2300 "logic error 1 in Step 9 of hmw_act");
2305 size_t n =
m_kk*j + k;
2307 sum2 += molality[j]*molality[k]*
m_PhiPhi_IJ[counterIJ];
2308 for (
size_t m = 1; m <
m_kk; m++) {
2312 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk[n];
2321 for (
size_t k = j+1; k <
m_kk; k++) {
2324 throw CanteraError(
"HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2325 "logic error 2 in Step 9 of hmw_act");
2330 size_t n =
m_kk*j + k;
2332 sum3 += molality[j]*molality[k]*
m_PhiPhi_IJ[counterIJ];
2333 for (
size_t m = 1; m <
m_kk; m++) {
2336 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk[n];
2345 for (
size_t k = 1; k <
m_kk; k++) {
2355 }
else if (k == j) {
2356 sum6 += 0.5 * molality[j]*molality[k]*
m_Lambda_nj(j,k);
2361 for (
size_t m = 1; m <
m_kk; m++) {
2367 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2373 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn[j];
2376 double sum_m_phi_minus_1 = 2.0 *
2377 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2380 double osmotic_coef;
2381 if (molalitysumUncropped > 1.0E-150) {
2382 osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2386 double lnwateract = -(
m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2413 for (
size_t k = 1; k <
m_kk; k++) {
2437 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2444 double molarcharge = 0.0;
2448 double molalitysum = 0.0;
2454 for (
size_t n = 1; n <
m_kk; n++) {
2458 molarcharge += fabs(
charge(n)) * molality[n];
2459 molalitysum += molality[n];
2473 for (
int z1 = 1; z1 <=4; z1++) {
2474 for (
int z2 =1; z2 <=4; z2++) {
2475 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
2482 for (
size_t i = 1; i < (
m_kk - 1); i++) {
2483 for (
size_t j = (i+1); j <
m_kk; j++) {
2485 size_t n =
m_kk*i + j;
2491 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
2492 if (x1 > 1.0E-100) {
2493 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2495 (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
2503 if (x2 > 1.0E-100) {
2504 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2506 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2522 for (
size_t i = 1; i <
m_kk - 1; i++) {
2523 for (
size_t j = i+1; j <
m_kk; j++) {
2525 size_t n =
m_kk*i + j;
2534 if (Is > 1.0E-150) {
2550 for (
size_t i = 1; i <
m_kk-1; i++) {
2551 for (
size_t j = i+1; j <
m_kk; j++) {
2553 size_t n =
m_kk*i + j;
2568 for (
size_t i = 1; i <
m_kk-1; i++) {
2569 for (
size_t j = i+1; j <
m_kk; j++) {
2571 size_t n =
m_kk*i + j;
2590 double dAphidT = dA_DebyedT /3.0;
2591 double dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2592 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2593 for (
size_t i = 1; i <
m_kk-1; i++) {
2594 for (
size_t j = i+1; j <
m_kk; j++) {
2596 size_t n =
m_kk*i + j;
2613 for (
size_t i = 1; i <
m_kk; i++) {
2623 for (
size_t j = 1; j <
m_kk; j++) {
2625 size_t n =
m_kk*i + j;
2630 sum1 += molality[j]*
2636 for (
size_t k = j+1; k <
m_kk; k++) {
2649 sum2 += molality[j]*(2.0*
m_Phi_IJ_L[counterIJ]);
2651 for (
size_t k = 1; k <
m_kk; k++) {
2661 molality[j]*molality[k]*
m_CMX_IJ_L[counterIJ2];
2672 for (
size_t k = 1; k <
m_kk; k++) {
2678 if (zeta_L != 0.0) {
2679 sum5 += molality[j]*molality[k]*zeta_L;
2688 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2701 for (
size_t j = 1; j <
m_kk; j++) {
2703 size_t n =
m_kk*i + j;
2708 sum1 += molality[j]*
2711 for (
size_t k = j+1; k <
m_kk; k++) {
2725 sum2 += molality[j]*(2.0*
m_Phi_IJ_L[counterIJ]);
2727 for (
size_t k = 1; k <
m_kk; k++) {
2735 sum4 += fabs(
charge(i)) *
2736 molality[j]*molality[k]*
m_CMX_IJ_L[counterIJ2];
2744 for (
size_t k = 1; k <
m_kk; k++) {
2751 if (zeta_L != 0.0) {
2752 sum5 += molality[j]*molality[k]*zeta_L;
2759 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2769 for (
size_t j = 1; j <
m_kk; j++) {
2773 for (
size_t k = 1; k <
m_kk; k++) {
2781 double sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_L[i];
2801 double term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
2803 for (
size_t j = 1; j <
m_kk; j++) {
2806 for (
size_t k = 1; k <
m_kk; k++) {
2809 size_t n =
m_kk*j + k;
2811 sum1 += molality[j]*molality[k]*
2816 for (
size_t k = j+1; k <
m_kk; k++) {
2817 if (j == (
m_kk-1)) {
2819 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2820 "logic error 1 in Step 9 of hmw_act");
2825 size_t n =
m_kk*j + k;
2828 for (
size_t m = 1; m <
m_kk; m++) {
2832 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_L[n];
2841 for (
size_t k = j+1; k <
m_kk; k++) {
2844 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2845 "logic error 2 in Step 9 of hmw_act");
2850 size_t n =
m_kk*j + k;
2853 for (
size_t m = 1; m <
m_kk; m++) {
2856 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_L[n];
2865 for (
size_t k = 1; k <
m_kk; k++) {
2875 }
else if (k == j) {
2881 for (
size_t m = 1; m <
m_kk; m++) {
2886 if (zeta_L != 0.0) {
2887 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
2893 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn_L[j];
2896 double sum_m_phi_minus_1 = 2.0 *
2897 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2900 double d_osmotic_coef_dT;
2901 if (molalitysum > 1.0E-150) {
2902 d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
2904 d_osmotic_coef_dT = 0.0;
2907 double d_lnwateract_dT = -(
m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
2932 for (
size_t k = 1; k <
m_kk; k++) {
2951 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2958 double molarcharge = 0.0;
2962 double molalitysum = 0.0;
2968 for (
size_t n = 1; n <
m_kk; n++) {
2972 molarcharge += fabs(
charge(n)) * molality[n];
2973 molalitysum += molality[n];
2987 for (
int z1 = 1; z1 <=4; z1++) {
2988 for (
int z2 =1; z2 <=4; z2++) {
2989 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
2996 for (
size_t i = 1; i < (
m_kk - 1); i++) {
2997 for (
size_t j = (i+1); j <
m_kk; j++) {
2999 size_t n =
m_kk*i + j;
3005 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
3006 if (x1 > 1.0E-100) {
3007 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
3009 (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
3017 if (x2 > 1.0E-100) {
3018 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3020 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3036 for (
size_t i = 1; i <
m_kk - 1; i++) {
3037 for (
size_t j = i+1; j <
m_kk; j++) {
3039 size_t n =
m_kk*i + j;
3048 if (Is > 1.0E-150) {
3064 for (
size_t i = 1; i <
m_kk-1; i++) {
3065 for (
size_t j = i+1; j <
m_kk; j++) {
3067 size_t n =
m_kk*i + j;
3082 for (
size_t i = 1; i <
m_kk-1; i++) {
3083 for (
size_t j = i+1; j <
m_kk; j++) {
3085 size_t n =
m_kk*i + j;
3104 double d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3105 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3106 for (
size_t i = 1; i <
m_kk-1; i++) {
3107 for (
size_t j = i+1; j <
m_kk; j++) {
3109 size_t n =
m_kk*i + j;
3121 d2FdT2 += molality[i]*molality[j] *
m_Phiprime_IJ[counterIJ];
3126 for (
size_t i = 1; i <
m_kk; i++) {
3136 for (
size_t j = 1; j <
m_kk; j++) {
3138 size_t n =
m_kk*i + j;
3143 sum1 += molality[j]*
3149 for (
size_t k = j+1; k <
m_kk; k++) {
3164 for (
size_t k = 1; k <
m_kk; k++) {
3173 sum4 += fabs(
charge(i)) *
3183 for (
size_t k = 1; k <
m_kk; k++) {
3189 if (zeta_LL != 0.0) {
3190 sum5 += molality[j]*molality[k]*zeta_LL;
3199 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3211 for (
size_t j = 1; j <
m_kk; j++) {
3213 size_t n =
m_kk*i + j;
3218 sum1 += molality[j]*
3221 for (
size_t k = j+1; k <
m_kk; k++) {
3237 for (
size_t k = 1; k <
m_kk; k++) {
3245 sum4 += fabs(
charge(i)) *
3255 for (
size_t k = 1; k <
m_kk; k++) {
3262 if (zeta_LL != 0.0) {
3263 sum5 += molality[j]*molality[k]*zeta_LL;
3270 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3279 for (
size_t j = 1; j <
m_kk; j++) {
3283 for (
size_t k = 1; k <
m_kk; k++) {
3291 double sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_LL[i];
3310 double term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3312 for (
size_t j = 1; j <
m_kk; j++) {
3315 for (
size_t k = 1; k <
m_kk; k++) {
3318 size_t n =
m_kk*j + k;
3321 sum1 += molality[j]*molality[k] *
3326 for (
size_t k = j+1; k <
m_kk; k++) {
3327 if (j == (
m_kk-1)) {
3329 throw CanteraError(
"HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3330 "logic error 1 in Step 9 of hmw_act");
3335 size_t n =
m_kk*j + k;
3338 for (
size_t m = 1; m <
m_kk; m++) {
3342 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_LL[n];
3351 for (
size_t k = j+1; k <
m_kk; k++) {
3354 throw CanteraError(
"HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3355 "logic error 2 in Step 9 of hmw_act");
3360 size_t n =
m_kk*j + k;
3364 for (
size_t m = 1; m <
m_kk; m++) {
3367 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_LL[n];
3376 for (
size_t k = 1; k <
m_kk; k++) {
3386 }
else if (k == j) {
3392 for (
size_t m = 1; m <
m_kk; m++) {
3397 if (zeta_LL != 0.0) {
3398 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
3405 sum7 += molality[j] * molality[j] * molality[j] *
m_Mu_nnn_LL[j];
3408 double sum_m_phi_minus_1 = 2.0 *
3409 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3412 double d2_osmotic_coef_dT2;
3413 if (molalitysum > 1.0E-150) {
3414 d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3416 d2_osmotic_coef_dT2 = 0.0;
3418 double d2_lnwateract_dT2 = -(
m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
3440 for (
size_t k = 1; k <
m_kk; k++) {
3458 double etheta[5][5], etheta_prime[5][5], sqrtIs;
3465 double molarcharge = 0.0;
3469 double molalitysum = 0.0;
3477 for (
size_t n = 1; n <
m_kk; n++) {
3481 molarcharge += fabs(
charge(n)) * molality[n];
3482 molalitysum += molality[n];
3496 for (
int z1 = 1; z1 <=4; z1++) {
3497 for (
int z2 =1; z2 <=4; z2++) {
3498 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
3505 for (
size_t i = 1; i < (
m_kk - 1); i++) {
3506 for (
size_t j = (i+1); j <
m_kk; j++) {
3508 size_t n =
m_kk*i + j;
3514 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
3515 if (x1 > 1.0E-100) {
3516 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
3518 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
3526 if (x2 > 1.0E-100) {
3527 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3529 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3545 for (
size_t i = 1; i <
m_kk - 1; i++) {
3546 for (
size_t j = i+1; j <
m_kk; j++) {
3548 size_t n =
m_kk*i + j;
3557 if (Is > 1.0E-150) {
3573 for (
size_t i = 1; i <
m_kk-1; i++) {
3574 for (
size_t j = i+1; j <
m_kk; j++) {
3576 size_t n =
m_kk*i + j;
3591 for (
size_t i = 1; i <
m_kk-1; i++) {
3592 for (
size_t j = i+1; j <
m_kk; j++) {
3594 size_t n =
m_kk*i + j;
3613 double dAphidP = dA_DebyedP /3.0;
3614 double dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3615 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3616 for (
size_t i = 1; i <
m_kk-1; i++) {
3617 for (
size_t j = i+1; j <
m_kk; j++) {
3619 size_t n =
m_kk*i + j;
3636 for (
size_t i = 1; i <
m_kk; i++) {
3646 for (
size_t j = 1; j <
m_kk; j++) {
3648 size_t n =
m_kk*i + j;
3653 sum1 += molality[j]*
3659 for (
size_t k = j+1; k <
m_kk; k++) {
3672 sum2 += molality[j]*(2.0*
m_Phi_IJ_P[counterIJ]);
3674 for (
size_t k = 1; k <
m_kk; k++) {
3683 sum4 += fabs(
charge(i)) *
3684 molality[j]*molality[k]*
m_CMX_IJ_P[counterIJ2];
3693 for (
size_t k = 1; k <
m_kk; k++) {
3699 if (zeta_P != 0.0) {
3700 sum5 += molality[j]*molality[k]*zeta_P;
3710 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3722 for (
size_t j = 1; j <
m_kk; j++) {
3724 size_t n =
m_kk*i + j;
3729 sum1 += molality[j] *
3732 for (
size_t k = j+1; k <
m_kk; k++) {
3746 sum2 += molality[j]*(2.0*
m_Phi_IJ_P[counterIJ]);
3748 for (
size_t k = 1; k <
m_kk; k++) {
3757 molality[j]*molality[k]*
m_CMX_IJ_P[counterIJ2];
3766 for (
size_t k = 1; k <
m_kk; k++) {
3773 if (zeta_P != 0.0) {
3774 sum5 += molality[j]*molality[k]*zeta_P;
3781 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3788 for (
size_t j = 1; j <
m_kk; j++) {
3792 for (
size_t k = 1; k <
m_kk; k++) {
3800 double sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_P[i];
3819 double term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3821 for (
size_t j = 1; j <
m_kk; j++) {
3824 for (
size_t k = 1; k <
m_kk; k++) {
3827 size_t n =
m_kk*j + k;
3829 sum1 += molality[j]*molality[k]*
3834 for (
size_t k = j+1; k <
m_kk; k++) {
3835 if (j == (
m_kk-1)) {
3837 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3838 "logic error 1 in Step 9 of hmw_act");
3843 size_t n =
m_kk*j + k;
3846 for (
size_t m = 1; m <
m_kk; m++) {
3850 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_P[n];
3859 for (
size_t k = j+1; k <
m_kk; k++) {
3862 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3863 "logic error 2 in Step 9 of hmw_act");
3868 size_t n =
m_kk*j + k;
3872 for (
size_t m = 1; m <
m_kk; m++) {
3875 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_P[n];
3884 for (
size_t k = 1; k <
m_kk; k++) {
3894 }
else if (k == j) {
3900 for (
size_t m = 1; m <
m_kk; m++) {
3905 if (zeta_P != 0.0) {
3906 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
3913 sum7 += molality[j] * molality[j] * molality[j] *
m_Mu_nnn_P[j];
3916 double sum_m_phi_minus_1 = 2.0 *
3917 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3921 double d_osmotic_coef_dP;
3922 if (molalitysum > 1.0E-150) {
3923 d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3925 d_osmotic_coef_dP = 0.0;
3927 double d_lnwateract_dP = -(
m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
3940 if( m_last_is == is ) {
3947 double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
3948 double aphi = 0.392;
3949 if (is < 1.0E-150) {
3950 for (
int i = 0; i < 17; i++) {
3959 for (
int i=1; i<=4; i++) {
3960 for (
int j=i; j<=4; j++) {
3964 double zprod = (double)ij;
3967 double x = 6.0* zprod * aphi * sqrt(is);
3969 double jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4)));
3971 double t = c3 * c4 * pow(x,c4);
3972 double dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
3973 double jprime = (jfunc/x)*(1.0 + jfunc*dj);
3975 elambda[ij] = zprod*jfunc / (4.0*is);
3976 elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
3983 double* etheta,
double* etheta_prime)
const 3990 "we shouldn't be here");
3992 "called with one species being neutral");
3998 *etheta_prime = 0.0;
4001 double f1 = (double)i / (2.0 * j);
4002 double f2 = (double)j / (2.0 * i);
4017 for (
size_t k = 1; k <
m_kk; k++) {
4024 double eterm = std::exp(-xoverc);
4026 double fptmp = IMS_bfCut_ - IMS_afCut_ /
IMS_cCut_ - IMS_bfCut_*xoverc
4027 + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
4028 double f_prime = 1.0 + eterm*fptmp;
4029 double f = xmolSolvent + IMS_efCut_
4030 + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
4032 double gptmp = IMS_bgCut_ - IMS_agCut_ /
IMS_cCut_ - IMS_bgCut_*xoverc
4033 + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
4034 double g_prime = 1.0 + eterm*gptmp;
4035 double g = xmolSolvent + IMS_egCut_
4036 + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
4038 double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
4039 double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
4040 double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
4042 tmp = log(xx) + lngammak;
4043 for (
size_t k = 1; k <
m_kk; k++) {
4060 writelog(
"Index Name MoleF MolalityCropped Charge\n");
4061 for (
size_t k = 0; k <
m_kk; k++) {
4062 writelogf(
"%2d %-16s %14.7le %14.7le %5.1f \n",
4066 writelog(
"\n Species Species beta0MX " 4067 "beta1MX beta2MX CphiMX alphaMX thetaij\n");
4068 for (
size_t i = 1; i <
m_kk - 1; i++) {
4069 for (
size_t j = i+1; j <
m_kk; j++) {
4070 size_t n = i *
m_kk + j;
4072 writelogf(
" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
4080 writelog(
"\n Species Species Species psi \n");
4081 for (
size_t i = 1; i <
m_kk; i++) {
4082 for (
size_t j = 1; j <
m_kk; j++) {
4083 for (
size_t k = 1; k <
m_kk; k++) {
4086 writelogf(
" %-16s %-16s %-16s %9.5f \n",
4103 doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4104 for (
size_t k = 0; k <
m_kk; k++) {
4105 acMolality[k] *= exp(
charge(k) * afac);
4118 doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4119 for (
size_t k = 0; k <
m_kk; k++) {
4133 doublereal afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
4134 for (
size_t k = 0; k <
m_kk; k++) {
4148 doublereal afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
4149 for (
size_t k = 0; k <
m_kk; k++) {
4163 doublereal afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
4164 for (
size_t k = 0; k <
m_kk; k++) {
4173 doublereal lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
4174 return lnGammaClMs2;
4181 return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
4188 return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
4195 return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
vector_fp m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P. Vector index is counterIJ.
vector_fp m_BMX_IJ_L
Derivative of BMX_IJ wrt T. Vector index is counterIJ.
vector_fp m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P. Vector index is counterIJ.
virtual void updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
void calcMolalitiesCropped() const
Calculate the cropped molalities.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
vector_fp m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT. Vector index is counterIJ.
doublereal molarVolume() const
Molar volume (m^3/kmol).
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
vector_fp m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T. Vector index is counterIJ.
void readXMLBinarySalt(XML_Node &BinSalt)
Process an XML node called "binarySaltParameters".
vector_fp m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P. Vector index is counterIJ.
vector_fp m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
int getId()
Get a unique id for a cached value.
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
vector_fp m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T. Vector index is counterIJ.
std::string name() const
Returns the name of the XML node.
vector_fp m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
virtual doublereal thermalExpansionCoeff() const
Return the volumetric thermal expansion coefficient. Units: 1/K.
doublereal MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
virtual void applyphScale(doublereal *acMolality) const
Apply the current phScale to a set of activity Coefficients or activities.
vector_fp IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
doublereal temperature() const
Temperature (K).
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
void readXMLLambdaNeutral(XML_Node &BinSalt)
Process an XML node called "lambdaNeutral".
size_t getFloatArray(const XML_Node &node, vector_fp &v, const bool convert, const std::string &unitsString, const std::string &nodeName)
This function reads the current node or a child node of the current node with the default name...
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
vector_fp m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT. Vector index is counterIJ.
doublereal m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
doublereal s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
vector_fp m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
vector_fp m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer's papers.
vector_fp m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P. Vector index is counterIJ.
vector_fp m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
const size_t npos
index returned by functions to indicate "no position"
void readXMLMunnnNeutral(XML_Node &BinSalt)
Process an XML node called "MunnnNeutral".
vector_fp m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer's paper.
vector_fp m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T. Vector index is counterIJ.
Header file for a common definitions used in electrolytes thermodynamics.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
vector_fp m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
virtual doublereal satPressure(doublereal T)
saturation pressure
virtual void setMolarDensity(const doublereal conc)
Set the internally stored molar density (kmol/m^3) for the phase.
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
size_t m_indexCLM
Index of the phScale species.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
vector_fp m_BMX_IJ_P
Derivative of BMX_IJ wrt P. Vector index is counterIJ.
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
virtual doublereal isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
vector_fp m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T. Vector index is counterIJ.
double elambda[17]
This is elambda, MEC.
Class XML_Node is a tree-based representation of the contents of an XML file.
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
vector_fp m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
vector_fp m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
doublereal s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale...
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the species standard states at the cu...
virtual void setDensity(const doublereal rho)
Set the internally stored density (kg/m^3) of the phase.
vector_fp m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P. Vector index is counterIJ.
void calcMCCutoffParams_()
Calculate molality cut-off parameters.
vector_fp m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P. Vector index is counterIJ.
virtual doublereal density() const
Density (kg/m^3).
void calcIMSCutoffParams_()
Precalculate the IMS Cutoff parameters for typeCutoff = 2.
vector_fp m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
void readXMLZetaCation(const XML_Node &BinSalt)
Process an XML node called "zetaCation".
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Initialize the phase parameters from an XML file.
doublereal m_xmolSolventMIN
vector_fp m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
vector_fp m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
vector_fp m_hfunc_IJ
hfunc, was called gprime in Pitzer's paper.
double elambda1[17]
This is elambda1, MEC.
vector_fp m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P. Vector index is counterIJ.
vector_fp m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT. Vector index is counterIJ.
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_V.
vector_fp m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT. Vector index is counterIJ.
doublereal IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
vector_fp m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
vector_fp m_Phi_IJ
Intermediate variable called Phi in Pitzer's paper.
vector_fp m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer's paper.
void setA_Debye(double A)
Set the A_Debye parameter.
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
int sign(T x)
Sign of a number. Returns -1 if x < 0, 1 if x > 0 and 0 if x == 0.
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure...
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
vector_fp m_Beta1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
vector_fp m_tmpV
vector of size m_kk, used as a temporary holding area.
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
HMWSoln()
Default Constructor.
vector_fp m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
vector_fp m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P. Vector index is counterIJ.
PDSS * m_waterSS
Water standard state calculator.
void readXMLTheta(XML_Node &BinSalt)
Process an XML node called "thetaAnion" or "thetaCation".
vector_fp m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer's paper.
vector_fp m_h2func_IJ
hfunc2, was called gprime in Pitzer's paper.
void printCoeffs() const
Print out all of the input Pitzer coefficients.
vector_fp m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
vector_fp m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T. Vector index is counterIJ.
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
vector_fp m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
const std::vector< XML_Node * > & children() const
Return an unchangeable reference to the vector of children of the current node.
doublereal s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
Headers for the HMWSoln ThermoPhase object, which models concentrated electrolyte solutions (see Ther...
std::unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
virtual doublereal relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
virtual void initThermo()
vector_fp m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
vector_fp m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
The WaterProps class is used to house several approximation routines for properties of water...
vector_fp m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT. Vector index is counterIJ.
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients...
CachedScalar getScalar(int id)
Get a reference to a CachedValue object representing a scalar (doublereal) with the given id...
vector_int CROP_speciesCropped_
This is a boolean-type vector indicating whether a species's activity coefficient is in the cropped r...
std::string speciesName(size_t k) const
Name of the species with index k.
vector_fp m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT. Vector index is counterIJ.
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer's papers.
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
doublereal IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
vector_fp m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT. Vector index is counterIJ.
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients' dependence on the temperature.
const int PHSCALE_NBS
Scale to be used for evaluation of single-ion activity coefficients is that used by the NBS standard ...
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale...
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
bool caseInsensitiveEquals(const std::string &input, const std::string &test)
Case insensitive equality predicate.
vector_int m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species...
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Base class for exceptions thrown by Cantera classes.
vector_fp m_CMX_IJ
Intermediate variable called CMX in Pitzer's paper.
int stateMFNumber() const
Return the State Mole Fraction Number.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
vector_fp m_CphiMX_ij_L
Derivative of Cphi_ij[i][j] wrt T. Vector index is counterIJ.
const int PHSCALE_PITZER
Scale to be used for the output of single-ion activity coefficients is that used by Pitzer...
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values.
vector_fp m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T. Vector index is counterIJ.
void readXMLPsi(XML_Node &BinSalt)
Process an XML node called "psiCommonAnion" or "psiCommonCation".
vector_fp m_g2func_IJ
This is the value of g2(x2) in Pitzer's papers. Vector index is counterIJ.
int m_pHScalingType
Scaling to be used for output of single-ion species activity coefficients.
virtual double d2A_DebyedT2_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the 2nd derivative of the Debye Huckel constant with respect to temperature as a function of...
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
void importPhase(XML_Node &phase, ThermoPhase *th)
Import a phase information into an empty ThermoPhase object.
vector_fp m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT. Vector index is counterIJ.
virtual doublereal pressure() const
Returns the current pressure of the phase.
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized activity concentrations.
vector_fp m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
vector_fp m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
vector_fp m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T. Vector index is counterIJ.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
XML_Node & child(const size_t n) const
Return a changeable reference to the n'th child of the current node.
#define AssertTrace(expr)
Assertion must be true or an error is thrown.
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
vector_fp m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT. Vector index is counterIJ.
vector_fp m_molalities
Current value of the molalities of the species in the phase.
double m_A_Debye
A_Debye: this expression appears on the top of the ln actCoeff term in the general Debye-Huckel expre...
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
void initLengths()
Initialize all of the species-dependent lengths in the object.
const doublereal SmallNumber
smallest number to compare to zero.
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
virtual doublereal relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.
virtual void initThermo()
void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
vector_fp m_gamma_tmp
Intermediate storage of the activity coefficient itself.
std::string id() const
Return the id attribute, if present.
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
Contains declarations for string manipulation functions within Cantera.
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
virtual void getActivities(doublereal *ac) const
Get the array of non-dimensional activities at the current solution temperature, pressure, and solution concentration.
doublereal m_weightSolvent
Molecular weight of the Solvent.
doublereal IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
vector_fp m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients...
vector_fp m_BMX_IJ
Intermediate variable called BMX in Pitzer's paper.
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature...
size_t m_kk
Number of species in the phase.
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
virtual double dA_DebyedP_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to pressure, as a function of tempe...
virtual double dA_DebyedT_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to temperature as a function of tem...
void s_updateScaling_pHScaling_dT2() const
Apply the current phScale to a set of 2nd derivatives of the activity Coefficients wrt temperature...
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
vector_fp m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer's paper.
T value
The value of the cached property.
vector_fp m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
Namespace for the Cantera kernel.
vector_fp m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model...
void calc_lambdas(double is) const
Calculate the lambda interactions.
vector_fp m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T. Vector index is counterIJ.
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_J.
vector_fp m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT. Vector index is counterIJ.
doublereal fpValueCheck(const std::string &val)
Translate a string into one doublereal value, with error checking.
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
virtual doublereal satPressure(doublereal T)
Get the saturation pressure for a given temperature.
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
vector_fp m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
vector_fp m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P. Vector index is counterIJ.
virtual double A_Debye_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the Debye Huckel constant as a function of temperature and pressure.
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of the species standard states at the current T and P of the solution...
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer's papers.
vector_fp m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T. Vector index is counterIJ.
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_L.
void setMoleFSolventMin(doublereal xmolSolventMIN)
Sets the minimum mole fraction in the molality formulation.
vector_fp m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P. Vector index is counterIJ.
vector_fp m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT. Vector index is counterIJ.
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
bool getOptionalFloat(const XML_Node &parent, const std::string &name, doublereal &fltRtn, const std::string &type)
Get an optional floating-point value from a child element.