28 double A_Debye_default = 1.172576;
29 double maxIionicStrength_default = 100.0;
30 double crop_ln_gamma_o_min_default = -6.0;
31 double crop_ln_gamma_o_max_default = 3.0;
32 double crop_ln_gamma_k_min_default = -5.0;
33 double crop_ln_gamma_k_max_default = 15.0;
42 m_maxIionicStrength(maxIionicStrength_default),
43 CROP_ln_gamma_o_min(crop_ln_gamma_o_min_default),
44 CROP_ln_gamma_o_max(crop_ln_gamma_o_max_default),
45 CROP_ln_gamma_k_min(crop_ln_gamma_k_min_default),
46 CROP_ln_gamma_k_max(crop_ln_gamma_k_max_default),
65 for (
size_t k = 0; k <
m_kk; k++) {
77 size_t kcation =
npos;
80 for (
size_t k = 0; k <
m_kk; k++) {
86 }
else if (
charge(k) < 0.0) {
93 if (kcation ==
npos || kanion ==
npos) {
96 double xuse = xcation;
98 if (xanion < xcation) {
100 if (
charge(kcation) != 1.0) {
104 if (
charge(kanion) != 1.0) {
108 xuse = xuse / factor;
137 return cp - beta * beta * tt * molarV / kappa_t;
166 for (
size_t k = 1; k <
m_kk; k++) {
175 double mvSolvent =
m_tmpV[0];
179 return 1.0 / mvSolvent;
191 for (
size_t k = 1; k <
m_kk; k++) {
204 for (
size_t k = 0; k <
m_kk; k++) {
205 acMolality[k] = exp(acMolality[k]);
223 for (
size_t k = 1; k <
m_kk; k++) {
237 for (
size_t k = 0; k <
m_kk; k++) {
245 for (
size_t k = 0; k <
m_kk; k++) {
257 for (
size_t k = 0; k <
m_kk; k++) {
268 for (
size_t k = 1; k <
m_kk; k++) {
280 for (
size_t k = 0; k <
m_kk; k++) {
293 for (
size_t k = 0; k <
m_kk; k++) {
301 for (
size_t k = 0; k <
m_kk; k++) {
310 for (
size_t k = 0; k <
m_kk; k++) {
328 static void check_nParams(
const string& method,
size_t nParams,
size_t m_formPitzerTemp)
330 if (m_formPitzerTemp == PITZER_TEMP_CONSTANT && nParams != 1) {
331 throw CanteraError(method,
"'constant' temperature model requires one"
332 " coefficient for each of parameter, but {} were given", nParams);
333 }
else if (m_formPitzerTemp == PITZER_TEMP_LINEAR && nParams != 2) {
334 throw CanteraError(method,
"'linear' temperature model requires two"
335 " coefficients for each parameter, but {} were given", nParams);
337 if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1 && nParams != 5) {
338 throw CanteraError(method,
"'complex' temperature model requires five"
339 " coefficients for each parameter, but {} were given", nParams);
343 void HMWSoln::setBinarySalt(
const string& sp1,
const string& sp2,
344 size_t nParams,
double* beta0,
double* beta1,
double* beta2,
345 double* Cphi,
double alpha1,
double alpha2)
350 throw CanteraError(
"HMWSoln::setBinarySalt",
"Species '{}' not found", sp1);
351 }
else if (k2 ==
npos) {
352 throw CanteraError(
"HMWSoln::setBinarySalt",
"Species '{}' not found", sp2);
357 throw CanteraError(
"HMWSoln::setBinarySalt",
"Species '{}' and '{}' "
358 "do not have opposite charges ({}, {})", sp1, sp2,
368 for (
size_t n = 0; n < nParams; n++) {
374 m_Alpha1MX_ij[c] = alpha1;
378 void HMWSoln::setTheta(
const string& sp1,
const string& sp2,
379 size_t nParams,
double* theta)
384 throw CanteraError(
"HMWSoln::setTheta",
"Species '{}' not found", sp1);
385 }
else if (k2 ==
npos) {
386 throw CanteraError(
"HMWSoln::setTheta",
"Species '{}' not found", sp2);
389 throw CanteraError(
"HMWSoln::setTheta",
"Species '{}' and '{}' "
390 "should both have the same (non-zero) charge ({}, {})", sp1, sp2,
396 for (
size_t n = 0; n < nParams; n++) {
401 void HMWSoln::setPsi(
const string& sp1,
const string& sp2,
402 const string& sp3,
size_t nParams,
double* psi)
408 throw CanteraError(
"HMWSoln::setPsi",
"Species '{}' not found", sp1);
409 }
else if (k2 ==
npos) {
410 throw CanteraError(
"HMWSoln::setPsi",
"Species '{}' not found", sp2);
411 }
else if (k3 ==
npos) {
412 throw CanteraError(
"HMWSoln::setPsi",
"Species '{}' not found", sp3);
417 throw CanteraError(
"HMWSoln::setPsi",
"All species must be ions and"
418 " must include at least one cation and one anion, but given species"
419 " (charges) were: {} ({}), {} ({}), and {} ({}).",
430 for (
size_t n = 0; n < nParams; n++) {
437 void HMWSoln::setLambda(
const string& sp1,
const string& sp2,
438 size_t nParams,
double* lambda)
443 throw CanteraError(
"HMWSoln::setLambda",
"Species '{}' not found", sp1);
444 }
else if (k2 ==
npos) {
445 throw CanteraError(
"HMWSoln::setLambda",
"Species '{}' not found", sp2);
449 throw CanteraError(
"HMWSoln::setLambda",
"Expected at least one neutral"
450 " species, but given species (charges) were: {} ({}) and {} ({}).",
457 size_t c = k1*
m_kk + k2;
458 for (
size_t n = 0; n < nParams; n++) {
464 void HMWSoln::setMunnn(
const string& sp,
size_t nParams,
double* munnn)
468 throw CanteraError(
"HMWSoln::setMunnn",
"Species '{}' not found", sp);
472 throw CanteraError(
"HMWSoln::setMunnn",
"Expected a neutral species,"
473 " got {} ({}).", sp,
charge(k));
476 for (
size_t n = 0; n < nParams; n++) {
482 void HMWSoln::setZeta(
const string& sp1,
const string& sp2,
483 const string& sp3,
size_t nParams,
double* psi)
489 throw CanteraError(
"HMWSoln::setZeta",
"Species '{}' not found", sp1);
490 }
else if (k2 ==
npos) {
491 throw CanteraError(
"HMWSoln::setZeta",
"Species '{}' not found", sp2);
492 }
else if (k3 ==
npos) {
493 throw CanteraError(
"HMWSoln::setZeta",
"Species '{}' not found", sp3);
498 throw CanteraError(
"HMWSoln::setZeta",
"Requires one neutral species, "
499 "one cation, and one anion, but given species (charges) were: "
500 "{} ({}), {} ({}), and {} ({}).",
507 }
else if (
charge(k3) == 0) {
519 for (
size_t n = 0; n < nParams; n++) {
525 void HMWSoln::setPitzerTempModel(
const string& model)
534 throw CanteraError(
"HMWSoln::setPitzerTempModel",
535 "Unknown Pitzer ActivityCoeff Temp model: {}", model);
549 void HMWSoln::setCroppingCoefficients(
double ln_gamma_k_min,
550 double ln_gamma_k_max,
double ln_gamma_o_min,
double ln_gamma_o_max)
552 CROP_ln_gamma_k_min = ln_gamma_k_min;
553 CROP_ln_gamma_k_max = ln_gamma_k_max;
554 CROP_ln_gamma_o_min = ln_gamma_o_min;
555 CROP_ln_gamma_o_max = ln_gamma_o_max;
558 vector<double> getSizedVector(
const AnyMap& item,
const string& key,
size_t nCoeffs)
561 if (item[key].is<double>()) {
564 v.push_back(item[key].asDouble());
566 v = item[key].asVector<
double>(1, nCoeffs);
568 if (v.size() == 1 && nCoeffs == 5) {
581 setPitzerTempModel(actData[
"temperature-model"].asString());
589 if (actData.hasKey(
"A_Debye")) {
590 if (actData[
"A_Debye"] ==
"variable") {
593 setA_Debye(actData.convert(
"A_Debye",
"kg^0.5/gmol^0.5"));
596 if (actData.hasKey(
"max-ionic-strength")) {
597 setMaxIonicStrength(actData[
"max-ionic-strength"].asDouble());
599 if (actData.hasKey(
"interactions")) {
600 for (
auto& item : actData[
"interactions"].asVector<AnyMap>()) {
601 auto&
species = item[
"species"].asVector<
string>(1, 3);
606 if (nsp == 2 && q0 * q1 < 0) {
608 vector<double> beta0 = getSizedVector(item,
"beta0", nCoeffs);
609 vector<double> beta1 = getSizedVector(item,
"beta1", nCoeffs);
610 vector<double> beta2 = getSizedVector(item,
"beta2", nCoeffs);
611 vector<double> Cphi = getSizedVector(item,
"Cphi", nCoeffs);
612 if (beta0.size() != beta1.size() || beta0.size() != beta2.size()
613 || beta0.size() != Cphi.size()) {
615 "Inconsistent binary salt array sizes ({}, {}, {}, {})",
616 beta0.size(), beta1.size(), beta2.size(), Cphi.size());
618 double alpha1 = item[
"alpha1"].asDouble();
619 double alpha2 = item.getDouble(
"alpha2", 0.0);
621 beta0.data(), beta1.data(), beta2.data(), Cphi.data(),
623 }
else if (nsp == 2 && q0 * q1 > 0) {
625 vector<double> theta = getSizedVector(item,
"theta", nCoeffs);
627 }
else if (nsp == 2 && q0 * q1 == 0) {
629 vector<double> lambda = getSizedVector(item,
"lambda", nCoeffs);
631 }
else if (nsp == 3 && q0 * q1 * q2 != 0) {
633 vector<double> psi = getSizedVector(item,
"psi", nCoeffs);
635 psi.size(), psi.data());
636 }
else if (nsp == 3 && q0 * q1 * q2 == 0) {
638 vector<double> zeta = getSizedVector(item,
"zeta", nCoeffs);
640 zeta.size(), zeta.data());
641 }
else if (nsp == 1) {
643 vector<double> mu = getSizedVector(item,
"mu", nCoeffs);
644 setMunnn(
species[0], mu.size(), mu.data());
648 if (actData.hasKey(
"cropping-coefficients")) {
649 auto& crop = actData[
"cropping-coefficients"].as<
AnyMap>();
650 setCroppingCoefficients(
651 crop.getDouble(
"ln_gamma_k_min", crop_ln_gamma_k_min_default),
652 crop.getDouble(
"ln_gamma_k_max", crop_ln_gamma_k_max_default),
653 crop.getDouble(
"ln_gamma_o_min", crop_ln_gamma_o_min_default),
654 crop.getDouble(
"ln_gamma_o_max", crop_ln_gamma_o_max_default));
660 for (
int i = 0; i < 17; i++) {
674 vector<double> mf(
m_kk, 0.0);
682 for (
size_t k = 0; k <
m_kk; k++) {
684 if (fabs(mf[k] *
charge(k)) > MaxC) {
691 if (fabs(sum) > 1.0E-30) {
693 if (mf[kHp] > sum * 1.1) {
707 if (mf[kOHm] > -sum * 1.1) {
719 if (notDone && kMaxC !=
npos) {
720 if (mf[kMaxC] > (1.1 * sum /
charge(kMaxC))) {
721 mf[kMaxC] -= sum /
charge(kMaxC);
722 mf[0] += sum /
charge(kMaxC);
741 void assignTrimmed(
AnyMap& interaction,
const string& key, vector<double>& values) {
742 while (values.size() > 1 && values.back() == 0) {
745 if (values.size() == 1) {
746 interaction[key] = values[0];
748 interaction[key] = values;
758 activityNode[
"temperature-model"] =
"linear";
761 activityNode[
"temperature-model"] =
"complex";
766 activityNode[
"A_Debye"] =
"variable";
767 }
else if (
m_A_Debye != A_Debye_default) {
768 activityNode[
"A_Debye"].setQuantity(
m_A_Debye,
"kg^0.5/gmol^0.5");
774 vector<AnyMap> interactions;
777 for (
size_t i = 1; i <
m_kk; i++) {
778 for (
size_t j = 1; j <
m_kk; j++) {
779 size_t c = i*
m_kk + j;
781 bool lambda_found =
false;
782 for (
size_t n = 0; n < nParams; n++) {
790 interaction[
"species"] = vector<string>{
792 vector<double> lambda(nParams);
793 for (
size_t n = 0; n < nParams; n++) {
796 assignTrimmed(interaction,
"lambda", lambda);
797 interactions.push_back(std::move(interaction));
802 if (c == 0 || i > j) {
807 bool salt_found =
false;
808 for (
size_t n = 0; n < nParams; n++) {
818 interaction[
"species"] = vector<string>{
820 vector<double> beta0(nParams), beta1(nParams), beta2(nParams), Cphi(nParams);
821 size_t last_nonzero = 0;
822 for (
size_t n = 0; n < nParams; n++) {
827 if (beta0[n] || beta1[n] || beta2[n] || Cphi[n]) {
831 if (last_nonzero == 0) {
832 interaction[
"beta0"] = beta0[0];
833 interaction[
"beta1"] = beta1[0];
834 interaction[
"beta2"] = beta2[0];
835 interaction[
"Cphi"] = Cphi[0];
837 beta0.resize(1 + last_nonzero);
838 beta1.resize(1 + last_nonzero);
839 beta2.resize(1 + last_nonzero);
840 Cphi.resize(1 + last_nonzero);
841 interaction[
"beta0"] = beta0;
842 interaction[
"beta1"] = beta1;
843 interaction[
"beta2"] = beta2;
844 interaction[
"Cphi"] = Cphi;
846 interaction[
"alpha1"] = m_Alpha1MX_ij[c];
850 interactions.push_back(std::move(interaction));
855 bool theta_found =
false;
856 for (
size_t n = 0; n < nParams; n++) {
864 interaction[
"species"] = vector<string>{
866 vector<double> theta(nParams);
867 for (
size_t n = 0; n < nParams; n++) {
870 assignTrimmed(interaction,
"theta", theta);
871 interactions.push_back(std::move(interaction));
880 for (
size_t i = 1; i <
m_kk; i++) {
884 for (
size_t j = i + 1; j <
m_kk; j++) {
888 for (
size_t k = j + 1; k <
m_kk; k++) {
893 for (
size_t n = 0; n < nParams; n++) {
896 interaction[
"species"] = vector<string>{
898 vector<double> psi(nParams);
899 for (
size_t m = 0; m < nParams; m++) {
902 assignTrimmed(interaction,
"psi", psi);
903 interactions.push_back(std::move(interaction));
912 for (
size_t i = 1; i <
m_kk; i++) {
916 for (
size_t j = 1; j <
m_kk; j++) {
920 for (
size_t k = 1; k <
m_kk; k++) {
925 for (
size_t n = 0; n < nParams; n++) {
928 interaction[
"species"] = vector<string>{
930 vector<double> zeta(nParams);
931 for (
size_t m = 0; m < nParams; m++) {
934 assignTrimmed(interaction,
"zeta", zeta);
935 interactions.push_back(std::move(interaction));
944 for (
size_t i = 1; i <
m_kk; i++) {
945 for (
size_t n = 0; n < nParams; n++) {
948 interaction[
"species"] = vector<string>{
speciesName(i)};
949 vector<double> mu(nParams);
950 for (
size_t m = 0; m < nParams; m++) {
953 assignTrimmed(interaction,
"mu", mu);
954 interactions.push_back(std::move(interaction));
960 activityNode[
"interactions"] = std::move(interactions);
963 if (CROP_ln_gamma_k_min != crop_ln_gamma_k_min_default) {
964 croppingCoeffs[
"ln_gamma_k_min"] = CROP_ln_gamma_k_min;
966 if (CROP_ln_gamma_k_max != crop_ln_gamma_k_max_default) {
967 croppingCoeffs[
"ln_gamma_k_max"] = CROP_ln_gamma_k_max;
969 if (CROP_ln_gamma_o_min != crop_ln_gamma_o_min_default) {
970 croppingCoeffs[
"ln_gamma_o_min"] = CROP_ln_gamma_o_min;
972 if (CROP_ln_gamma_o_max != crop_ln_gamma_o_max_default) {
973 croppingCoeffs[
"ln_gamma_o_max"] = CROP_ln_gamma_o_max;
975 if (croppingCoeffs.
size()) {
976 activityNode[
"cropping-coefficients"] = std::move(croppingCoeffs);
979 phaseNode[
"activity-data"] = std::move(activityNode);
986 if (tempArg != -1.0) {
990 if (presArg != -1.0) {
1009 throw CanteraError(
"HMWSoln::A_Debye_TP",
"shouldn't be here");
1017 if (tempArg != -1.0) {
1021 if (presArg != -1.0) {
1033 throw CanteraError(
"HMWSoln::dA_DebyedT_TP",
"shouldn't be here");
1041 if (tempArg != -1.0) {
1045 if (presArg != -1.0) {
1058 dAdP = cached.
value;
1061 cached.
value = dAdP;
1065 throw CanteraError(
"HMWSoln::dA_DebyedP_TP",
"shouldn't be here");
1073 double dAphidT = dAdT /3.0;
1075 if (tempArg != -1.0) {
1084 double dAphidP = dAdP /3.0;
1086 if (tempArg != -1.0) {
1095 if (tempArg != -1.0) {
1100 double d2Aphi = d2 / 3.0;
1101 return 2.0 * A_L / T + 4.0 *
GasConstant * T * T *d2Aphi;
1107 if (tempArg != -1.0) {
1111 if (presArg != -1.0) {
1123 throw CanteraError(
"HMWSoln::d2A_DebyedT2_TP",
"shouldn't be here");
1137 size_t maxCounterIJlen = 1 + (
m_kk-1) * (
m_kk-2) / 2;
1140 int TCoeffLength = 1;
1171 m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
1212 m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1224 m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1233 m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1273 double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1274 double lnxs = log(xx);
1276 for (
size_t k = 1; k <
m_kk; k++) {
1314 for (
size_t k = 0; k <
m_kk; k++) {
1320 if (cropMethod == 0) {
1327 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1328 double charge_i =
charge(i);
1329 double abs_charge_i = fabs(charge_i);
1330 if (charge_i == 0.0) {
1333 for (
size_t j = (i+1); j <
m_kk; j++) {
1334 double charge_j =
charge(j);
1335 double abs_charge_j = fabs(charge_j);
1338 if (charge_i * charge_j < 0) {
1343 if (Imax > Iac_max) {
1347 if (Imax > Iac_max) {
1352 if (Imax > Iac_max) {
1356 if (Imax > Iac_max) {
1366 for (
int times = 0; times< 10; times++) {
1367 double anion_charge = 0.0;
1368 double cation_charge = 0.0;
1369 size_t anion_contrib_max_i =
npos;
1370 double anion_contrib_max = -1.0;
1371 size_t cation_contrib_max_i =
npos;
1372 double cation_contrib_max = -1.0;
1373 for (
size_t i = 0; i <
m_kk; i++) {
1374 double charge_i =
charge(i);
1375 if (charge_i < 0.0) {
1377 anion_charge += anion_contrib;
1378 if (anion_contrib > anion_contrib_max) {
1379 anion_contrib_max = anion_contrib;
1380 anion_contrib_max_i = i;
1382 }
else if (charge_i > 0.0) {
1384 cation_charge += cation_contrib;
1385 if (cation_contrib > cation_contrib_max) {
1386 cation_contrib_max = cation_contrib;
1387 cation_contrib_max_i = i;
1391 double total_charge = cation_charge - anion_charge;
1392 if (total_charge > 1.0E-8) {
1393 double desiredCrop = total_charge/
charge(cation_contrib_max_i);
1395 if (desiredCrop < maxCrop) {
1401 }
else if (total_charge < -1.0E-8) {
1402 double desiredCrop = total_charge/
charge(anion_contrib_max_i);
1404 if (desiredCrop < maxCrop) {
1416 if (cropMethod == 1) {
1419 double xmolSolvent = molF[0];
1425 double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1426 double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1428 for (
size_t k = 0; k <
m_kk; k++) {
1436 for (
size_t k = 0; k <
m_kk; k++) {
1441 for (
size_t k = 0; k <
m_kk; k++) {
1454 for (
size_t i = 0; i <
m_kk; i++) {
1456 size_t nc =
m_kk * i;
1460 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1461 size_t n =
m_kk * i + i;
1463 for (
size_t j = (i+1); j <
m_kk; j++) {
1465 size_t nc =
m_kk * i + j;
1475 double IMS_gamma_o_min_ = 1.0E-5;
1476 double IMS_gamma_k_min_ = 10.0;
1477 double IMS_slopefCut_ = 0.6;
1479 IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
1481 bool converged =
false;
1483 for (
int its = 0; its < 100 && !converged; its++) {
1485 IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) -IMS_efCut_;
1486 IMS_bfCut_ = IMS_afCut_ /
IMS_cCut_ + IMS_slopefCut_ - 1.0;
1492 IMS_efCut_ = - eterm * tmp;
1493 if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
1499 " failed to converge on the f polynomial");
1502 double f_0 = IMS_afCut_ + IMS_efCut_;
1503 double f_prime_0 = 1.0 - IMS_afCut_ /
IMS_cCut_ + IMS_bfCut_;
1505 for (
int its = 0; its < 100 && !converged; its++) {
1507 double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
1508 IMS_agCut_ = exp(lng_0) - IMS_egCut_;
1515 IMS_egCut_ = - eterm * tmp;
1516 if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
1522 " failed to converge on the g polynomial");
1528 double MC_X_o_min_ = 0.35;
1530 double MC_slopepCut_ = 0.02;
1534 MC_apCut_ = MC_X_o_min_;
1536 bool converged =
false;
1539 for (
int its = 0; its < 500 && !converged; its++) {
1541 MC_apCut_ = damp *(MC_X_o_min_ - MC_epCut_) + (1-damp) * MC_apCut_;
1542 double MC_bpCutNew = MC_apCut_ / MC_cpCut_ + MC_slopepCut_ - 1.0;
1543 MC_bpCut_ = damp * MC_bpCutNew + (1-damp) * MC_bpCut_;
1544 double MC_dpCutNew = ((- MC_apCut_/MC_cpCut_ + MC_bpCut_ - MC_bpCut_ *
MC_X_o_cutoff_/MC_cpCut_)
1547 MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
1550 MC_epCut_ = - eterm * tmp;
1551 double diff = MC_epCut_ - oldV;
1552 if (fabs(diff) < 1.0E-14) {
1558 " failed to converge on the p polynomial");
1565 const double twoT = 2.0 * T;
1566 const double invT = 1.0 / T;
1567 const double invT2 = invT * invT;
1568 const double twoinvT3 = 2.0 * invT * invT2;
1569 double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
1579 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1580 for (
size_t j = (i+1); j <
m_kk; j++) {
1582 size_t n =
m_kk*i + j;
1592 case PITZER_TEMP_CONSTANT:
1594 case PITZER_TEMP_LINEAR:
1597 + beta0MX_coeff[1]*tlin;
1601 + beta1MX_coeff[1]*tlin;
1605 + beta2MX_coeff[1]*tlin;
1609 + CphiMX_coeff[1]*tlin;
1612 m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1617 case PITZER_TEMP_COMPLEX1:
1619 + beta0MX_coeff[1]*tlin
1620 + beta0MX_coeff[2]*tquad
1621 + beta0MX_coeff[3]*tinv
1622 + beta0MX_coeff[4]*tln;
1624 + beta1MX_coeff[1]*tlin
1625 + beta1MX_coeff[2]*tquad
1626 + beta1MX_coeff[3]*tinv
1627 + beta1MX_coeff[4]*tln;
1629 + beta2MX_coeff[1]*tlin
1630 + beta2MX_coeff[2]*tquad
1631 + beta2MX_coeff[3]*tinv
1632 + beta2MX_coeff[4]*tln;
1634 + CphiMX_coeff[1]*tlin
1635 + CphiMX_coeff[2]*tquad
1636 + CphiMX_coeff[3]*tinv
1637 + CphiMX_coeff[4]*tln;
1639 + Theta_coeff[1]*tlin
1640 + Theta_coeff[2]*tquad
1641 + Theta_coeff[3]*tinv
1642 + Theta_coeff[4]*tln;
1644 + beta0MX_coeff[2]*twoT
1645 - beta0MX_coeff[3]*invT2
1646 + beta0MX_coeff[4]*invT;
1648 + beta1MX_coeff[2]*twoT
1649 - beta1MX_coeff[3]*invT2
1650 + beta1MX_coeff[4]*invT;
1652 + beta2MX_coeff[2]*twoT
1653 - beta2MX_coeff[3]*invT2
1654 + beta2MX_coeff[4]*invT;
1656 + CphiMX_coeff[2]*twoT
1657 - CphiMX_coeff[3]*invT2
1658 + CphiMX_coeff[4]*invT;
1660 + Theta_coeff[2]*twoT
1661 - Theta_coeff[3]*invT2
1662 + Theta_coeff[4]*invT;
1666 + beta0MX_coeff[2]*2.0
1667 + beta0MX_coeff[3]*twoinvT3
1668 - beta0MX_coeff[4]*invT2;
1670 + beta1MX_coeff[2]*2.0
1671 + beta1MX_coeff[3]*twoinvT3
1672 - beta1MX_coeff[4]*invT2;
1674 + beta2MX_coeff[2]*2.0
1675 + beta2MX_coeff[3]*twoinvT3
1676 - beta2MX_coeff[4]*invT2;
1678 + CphiMX_coeff[2]*2.0
1679 + CphiMX_coeff[3]*twoinvT3
1680 - CphiMX_coeff[4]*invT2;
1682 + Theta_coeff[2]*2.0
1683 + Theta_coeff[3]*twoinvT3
1684 - Theta_coeff[4]*invT2;
1694 for (
size_t i = 1; i <
m_kk; i++) {
1696 for (
size_t j = 1; j <
m_kk; j++) {
1697 size_t n = i *
m_kk + j;
1700 case PITZER_TEMP_CONSTANT:
1703 case PITZER_TEMP_LINEAR:
1704 m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1708 case PITZER_TEMP_COMPLEX1:
1710 + Lambda_coeff[1]*tlin
1711 + Lambda_coeff[2]*tquad
1712 + Lambda_coeff[3]*tinv
1713 + Lambda_coeff[4]*tln;
1716 + Lambda_coeff[2]*twoT
1717 - Lambda_coeff[3]*invT2
1718 + Lambda_coeff[4]*invT;
1722 + Lambda_coeff[3]*twoinvT3
1723 - Lambda_coeff[4]*invT2;
1729 case PITZER_TEMP_CONSTANT:
1732 case PITZER_TEMP_LINEAR:
1733 m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1737 case PITZER_TEMP_COMPLEX1:
1749 + Mu_coeff[3]*twoinvT3
1750 - Mu_coeff[4]*invT2;
1758 case PITZER_TEMP_CONSTANT:
1759 for (
size_t i = 1; i <
m_kk; i++) {
1760 for (
size_t j = 1; j <
m_kk; j++) {
1761 for (
size_t k = 1; k <
m_kk; k++) {
1769 case PITZER_TEMP_LINEAR:
1770 for (
size_t i = 1; i <
m_kk; i++) {
1771 for (
size_t j = 1; j <
m_kk; j++) {
1772 for (
size_t k = 1; k <
m_kk; k++) {
1775 m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1782 case PITZER_TEMP_COMPLEX1:
1783 for (
size_t i = 1; i <
m_kk; i++) {
1784 for (
size_t j = 1; j <
m_kk; j++) {
1785 for (
size_t k = 1; k <
m_kk; k++) {
1790 + Psi_coeff[2]*tquad
1795 - Psi_coeff[3]*invT2
1796 + Psi_coeff[4]*invT;
1799 + Psi_coeff[3]*twoinvT3
1800 - Psi_coeff[4]*invT2;
1818 double etheta[5][5], etheta_prime[5][5], sqrtIs;
1825 double molarcharge = 0.0;
1829 double molalitysumUncropped = 0.0;
1835 for (
size_t n = 1; n <
m_kk; n++) {
1839 molarcharge += fabs(
charge(n)) * molality[n];
1854 for (
int z1 = 1; z1 <=4; z1++) {
1855 for (
int z2 =1; z2 <=4; z2++) {
1856 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
1863 for (
size_t i = 1; i < (
m_kk - 1); i++) {
1864 for (
size_t j = (i+1); j <
m_kk; j++) {
1866 size_t n =
m_kk*i + j;
1872 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
1873 if (x1 > 1.0E-100) {
1874 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
1876 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
1884 if (x2 > 1.0E-100) {
1885 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
1887 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
1902 for (
size_t i = 1; i <
m_kk - 1; i++) {
1903 for (
size_t j = i+1; j <
m_kk; j++) {
1905 size_t n =
m_kk*i + j;
1915 if (Is > 1.0E-150) {
1932 for (
size_t i = 1; i <
m_kk-1; i++) {
1933 for (
size_t j = i+1; j <
m_kk; j++) {
1935 size_t n =
m_kk*i + j;
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 int z1 = (int) fabs(
charge(i));
1961 int z2 = (int) fabs(
charge(j));
1976 double F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
1977 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
1978 for (
size_t i = 1; i <
m_kk-1; i++) {
1979 for (
size_t j = i+1; j <
m_kk; j++) {
1981 size_t n =
m_kk*i + j;
1998 for (
size_t i = 1; i <
m_kk; i++) {
2011 for (
size_t j = 1; j <
m_kk; j++) {
2013 size_t n =
m_kk*i + j;
2018 sum1 += molality[j] *
2024 for (
size_t k = j+1; k <
m_kk; k++) {
2028 sum3 += molality[j]*molality[k]*
m_Psi_ijk[n];
2037 sum2 += molality[j]*(2.0*
m_Phi_IJ[counterIJ]);
2039 for (
size_t k = 1; k <
m_kk; k++) {
2043 sum2 += molality[j]*molality[k]*
m_Psi_ijk[n];
2048 sum4 += (fabs(
charge(i))*
2049 molality[j]*molality[k]*
m_CMX_IJ[counterIJ2]);
2059 for (
size_t k = 1; k <
m_kk; k++) {
2066 sum5 += molality[j]*molality[k]*zeta;
2090 for (
size_t j = 1; j <
m_kk; j++) {
2092 size_t n =
m_kk*i + j;
2097 sum1 += molality[j]*
2100 for (
size_t k = j+1; k <
m_kk; k++) {
2104 sum3 += molality[j]*molality[k]*
m_Psi_ijk[n];
2114 sum2 += molality[j]*(2.0*
m_Phi_IJ[counterIJ]);
2116 for (
size_t k = 1; k <
m_kk; k++) {
2120 sum2 += molality[j]*molality[k]*
m_Psi_ijk[n];
2125 molality[j]*molality[k]*
m_CMX_IJ[counterIJ2];
2134 for (
size_t k = 1; k <
m_kk; k++) {
2142 sum5 += molality[j]*molality[k]*zeta;
2158 for (
size_t j = 1; j <
m_kk; j++) {
2162 for (
size_t k = 1; k <
m_kk; k++) {
2165 sum3 += molality[j]*molality[k]*
m_Psi_ijk[n];
2170 double sum2 = 3.0 * molality[i]* molality[i] *
m_Mu_nnn[i];
2192 double term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2194 for (
size_t j = 1; j <
m_kk; j++) {
2197 for (
size_t k = 1; k <
m_kk; k++) {
2200 size_t n =
m_kk*j + k;
2203 sum1 += molality[j]*molality[k]*
2208 for (
size_t k = j+1; k <
m_kk; k++) {
2209 if (j == (
m_kk-1)) {
2211 throw CanteraError(
"HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2212 "logic error 1 in Step 9 of hmw_act");
2217 size_t n =
m_kk*j + k;
2219 sum2 += molality[j]*molality[k]*
m_PhiPhi_IJ[counterIJ];
2220 for (
size_t m = 1; m <
m_kk; m++) {
2224 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk[n];
2233 for (
size_t k = j+1; k <
m_kk; k++) {
2236 throw CanteraError(
"HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2237 "logic error 2 in Step 9 of hmw_act");
2242 size_t n =
m_kk*j + k;
2244 sum3 += molality[j]*molality[k]*
m_PhiPhi_IJ[counterIJ];
2245 for (
size_t m = 1; m <
m_kk; m++) {
2248 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk[n];
2257 for (
size_t k = 1; k <
m_kk; k++) {
2267 }
else if (k == j) {
2268 sum6 += 0.5 * molality[j]*molality[k]*
m_Lambda_nj(j,k);
2273 for (
size_t m = 1; m <
m_kk; m++) {
2279 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2285 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn[j];
2288 double sum_m_phi_minus_1 = 2.0 *
2289 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2292 double osmotic_coef;
2293 if (molalitysumUncropped > 1.0E-150) {
2294 osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2298 double lnwateract = -(
m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2325 for (
size_t k = 1; k <
m_kk; k++) {
2349 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2356 double molarcharge = 0.0;
2360 double molalitysum = 0.0;
2366 for (
size_t n = 1; n <
m_kk; n++) {
2370 molarcharge += fabs(
charge(n)) * molality[n];
2371 molalitysum += molality[n];
2385 for (
int z1 = 1; z1 <=4; z1++) {
2386 for (
int z2 =1; z2 <=4; z2++) {
2387 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
2394 for (
size_t i = 1; i < (
m_kk - 1); i++) {
2395 for (
size_t j = (i+1); j <
m_kk; j++) {
2397 size_t n =
m_kk*i + j;
2403 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
2404 if (x1 > 1.0E-100) {
2405 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2407 (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
2415 if (x2 > 1.0E-100) {
2416 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2418 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2434 for (
size_t i = 1; i <
m_kk - 1; i++) {
2435 for (
size_t j = i+1; j <
m_kk; j++) {
2437 size_t n =
m_kk*i + j;
2446 if (Is > 1.0E-150) {
2462 for (
size_t i = 1; i <
m_kk-1; i++) {
2463 for (
size_t j = i+1; j <
m_kk; j++) {
2465 size_t n =
m_kk*i + j;
2480 for (
size_t i = 1; i <
m_kk-1; i++) {
2481 for (
size_t j = i+1; j <
m_kk; j++) {
2483 size_t n =
m_kk*i + j;
2502 double dAphidT = dA_DebyedT /3.0;
2503 double dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2504 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2505 for (
size_t i = 1; i <
m_kk-1; i++) {
2506 for (
size_t j = i+1; j <
m_kk; j++) {
2508 size_t n =
m_kk*i + j;
2525 for (
size_t i = 1; i <
m_kk; i++) {
2535 for (
size_t j = 1; j <
m_kk; j++) {
2537 size_t n =
m_kk*i + j;
2542 sum1 += molality[j]*
2548 for (
size_t k = j+1; k <
m_kk; k++) {
2561 sum2 += molality[j]*(2.0*
m_Phi_IJ_L[counterIJ]);
2563 for (
size_t k = 1; k <
m_kk; k++) {
2573 molality[j]*molality[k]*
m_CMX_IJ_L[counterIJ2];
2584 for (
size_t k = 1; k <
m_kk; k++) {
2590 if (zeta_L != 0.0) {
2591 sum5 += molality[j]*molality[k]*zeta_L;
2600 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2613 for (
size_t j = 1; j <
m_kk; j++) {
2615 size_t n =
m_kk*i + j;
2620 sum1 += molality[j]*
2623 for (
size_t k = j+1; k <
m_kk; k++) {
2637 sum2 += molality[j]*(2.0*
m_Phi_IJ_L[counterIJ]);
2639 for (
size_t k = 1; k <
m_kk; k++) {
2647 sum4 += fabs(
charge(i)) *
2648 molality[j]*molality[k]*
m_CMX_IJ_L[counterIJ2];
2656 for (
size_t k = 1; k <
m_kk; k++) {
2663 if (zeta_L != 0.0) {
2664 sum5 += molality[j]*molality[k]*zeta_L;
2671 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2681 for (
size_t j = 1; j <
m_kk; j++) {
2685 for (
size_t k = 1; k <
m_kk; k++) {
2693 double sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_L[i];
2713 double term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
2715 for (
size_t j = 1; j <
m_kk; j++) {
2718 for (
size_t k = 1; k <
m_kk; k++) {
2721 size_t n =
m_kk*j + k;
2723 sum1 += molality[j]*molality[k]*
2728 for (
size_t k = j+1; k <
m_kk; k++) {
2729 if (j == (
m_kk-1)) {
2731 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2732 "logic error 1 in Step 9 of hmw_act");
2737 size_t n =
m_kk*j + k;
2740 for (
size_t m = 1; m <
m_kk; m++) {
2744 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_L[n];
2753 for (
size_t k = j+1; k <
m_kk; k++) {
2756 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2757 "logic error 2 in Step 9 of hmw_act");
2762 size_t n =
m_kk*j + k;
2765 for (
size_t m = 1; m <
m_kk; m++) {
2768 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_L[n];
2777 for (
size_t k = 1; k <
m_kk; k++) {
2787 }
else if (k == j) {
2793 for (
size_t m = 1; m <
m_kk; m++) {
2798 if (zeta_L != 0.0) {
2799 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
2805 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn_L[j];
2808 double sum_m_phi_minus_1 = 2.0 *
2809 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2812 double d_osmotic_coef_dT;
2813 if (molalitysum > 1.0E-150) {
2814 d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
2816 d_osmotic_coef_dT = 0.0;
2819 double d_lnwateract_dT = -(
m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
2844 for (
size_t k = 1; k <
m_kk; k++) {
2863 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2870 double molarcharge = 0.0;
2874 double molalitysum = 0.0;
2880 for (
size_t n = 1; n <
m_kk; n++) {
2884 molarcharge += fabs(
charge(n)) * molality[n];
2885 molalitysum += molality[n];
2899 for (
int z1 = 1; z1 <=4; z1++) {
2900 for (
int z2 =1; z2 <=4; z2++) {
2901 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
2908 for (
size_t i = 1; i < (
m_kk - 1); i++) {
2909 for (
size_t j = (i+1); j <
m_kk; j++) {
2911 size_t n =
m_kk*i + j;
2917 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
2918 if (x1 > 1.0E-100) {
2919 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
2921 (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
2929 if (x2 > 1.0E-100) {
2930 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2932 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2948 for (
size_t i = 1; i <
m_kk - 1; i++) {
2949 for (
size_t j = i+1; j <
m_kk; j++) {
2951 size_t n =
m_kk*i + j;
2960 if (Is > 1.0E-150) {
2976 for (
size_t i = 1; i <
m_kk-1; i++) {
2977 for (
size_t j = i+1; j <
m_kk; j++) {
2979 size_t n =
m_kk*i + j;
2994 for (
size_t i = 1; i <
m_kk-1; i++) {
2995 for (
size_t j = i+1; j <
m_kk; j++) {
2997 size_t n =
m_kk*i + j;
3016 double d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3017 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3018 for (
size_t i = 1; i <
m_kk-1; i++) {
3019 for (
size_t j = i+1; j <
m_kk; j++) {
3021 size_t n =
m_kk*i + j;
3033 d2FdT2 += molality[i]*molality[j] *
m_Phiprime_IJ[counterIJ];
3038 for (
size_t i = 1; i <
m_kk; i++) {
3048 for (
size_t j = 1; j <
m_kk; j++) {
3050 size_t n =
m_kk*i + j;
3055 sum1 += molality[j]*
3061 for (
size_t k = j+1; k <
m_kk; k++) {
3076 for (
size_t k = 1; k <
m_kk; k++) {
3085 sum4 += fabs(
charge(i)) *
3095 for (
size_t k = 1; k <
m_kk; k++) {
3101 if (zeta_LL != 0.0) {
3102 sum5 += molality[j]*molality[k]*zeta_LL;
3111 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3123 for (
size_t j = 1; j <
m_kk; j++) {
3125 size_t n =
m_kk*i + j;
3130 sum1 += molality[j]*
3133 for (
size_t k = j+1; k <
m_kk; k++) {
3149 for (
size_t k = 1; k <
m_kk; k++) {
3157 sum4 += fabs(
charge(i)) *
3167 for (
size_t k = 1; k <
m_kk; k++) {
3174 if (zeta_LL != 0.0) {
3175 sum5 += molality[j]*molality[k]*zeta_LL;
3182 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3191 for (
size_t j = 1; j <
m_kk; j++) {
3195 for (
size_t k = 1; k <
m_kk; k++) {
3203 double sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_LL[i];
3222 double term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3224 for (
size_t j = 1; j <
m_kk; j++) {
3227 for (
size_t k = 1; k <
m_kk; k++) {
3230 size_t n =
m_kk*j + k;
3233 sum1 += molality[j]*molality[k] *
3238 for (
size_t k = j+1; k <
m_kk; k++) {
3239 if (j == (
m_kk-1)) {
3241 throw CanteraError(
"HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3242 "logic error 1 in Step 9 of hmw_act");
3247 size_t n =
m_kk*j + k;
3250 for (
size_t m = 1; m <
m_kk; m++) {
3254 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_LL[n];
3263 for (
size_t k = j+1; k <
m_kk; k++) {
3266 throw CanteraError(
"HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3267 "logic error 2 in Step 9 of hmw_act");
3272 size_t n =
m_kk*j + k;
3276 for (
size_t m = 1; m <
m_kk; m++) {
3279 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_LL[n];
3288 for (
size_t k = 1; k <
m_kk; k++) {
3298 }
else if (k == j) {
3304 for (
size_t m = 1; m <
m_kk; m++) {
3309 if (zeta_LL != 0.0) {
3310 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
3317 sum7 += molality[j] * molality[j] * molality[j] *
m_Mu_nnn_LL[j];
3320 double sum_m_phi_minus_1 = 2.0 *
3321 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3324 double d2_osmotic_coef_dT2;
3325 if (molalitysum > 1.0E-150) {
3326 d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3328 d2_osmotic_coef_dT2 = 0.0;
3330 double d2_lnwateract_dT2 = -(
m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
3352 for (
size_t k = 1; k <
m_kk; k++) {
3370 double etheta[5][5], etheta_prime[5][5], sqrtIs;
3377 double molarcharge = 0.0;
3381 double molalitysum = 0.0;
3389 for (
size_t n = 1; n <
m_kk; n++) {
3393 molarcharge += fabs(
charge(n)) * molality[n];
3394 molalitysum += molality[n];
3408 for (
int z1 = 1; z1 <=4; z1++) {
3409 for (
int z2 =1; z2 <=4; z2++) {
3410 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
3417 for (
size_t i = 1; i < (
m_kk - 1); i++) {
3418 for (
size_t j = (i+1); j <
m_kk; j++) {
3420 size_t n =
m_kk*i + j;
3426 double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
3427 if (x1 > 1.0E-100) {
3428 m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
3430 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
3438 if (x2 > 1.0E-100) {
3439 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3441 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3457 for (
size_t i = 1; i <
m_kk - 1; i++) {
3458 for (
size_t j = i+1; j <
m_kk; j++) {
3460 size_t n =
m_kk*i + j;
3469 if (Is > 1.0E-150) {
3485 for (
size_t i = 1; i <
m_kk-1; i++) {
3486 for (
size_t j = i+1; j <
m_kk; j++) {
3488 size_t n =
m_kk*i + j;
3503 for (
size_t i = 1; i <
m_kk-1; i++) {
3504 for (
size_t j = i+1; j <
m_kk; j++) {
3506 size_t n =
m_kk*i + j;
3525 double dAphidP = dA_DebyedP /3.0;
3526 double dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3527 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3528 for (
size_t i = 1; i <
m_kk-1; i++) {
3529 for (
size_t j = i+1; j <
m_kk; j++) {
3531 size_t n =
m_kk*i + j;
3548 for (
size_t i = 1; i <
m_kk; i++) {
3558 for (
size_t j = 1; j <
m_kk; j++) {
3560 size_t n =
m_kk*i + j;
3565 sum1 += molality[j]*
3571 for (
size_t k = j+1; k <
m_kk; k++) {
3584 sum2 += molality[j]*(2.0*
m_Phi_IJ_P[counterIJ]);
3586 for (
size_t k = 1; k <
m_kk; k++) {
3595 sum4 += fabs(
charge(i)) *
3596 molality[j]*molality[k]*
m_CMX_IJ_P[counterIJ2];
3605 for (
size_t k = 1; k <
m_kk; k++) {
3611 if (zeta_P != 0.0) {
3612 sum5 += molality[j]*molality[k]*zeta_P;
3622 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3634 for (
size_t j = 1; j <
m_kk; j++) {
3636 size_t n =
m_kk*i + j;
3641 sum1 += molality[j] *
3644 for (
size_t k = j+1; k <
m_kk; k++) {
3658 sum2 += molality[j]*(2.0*
m_Phi_IJ_P[counterIJ]);
3660 for (
size_t k = 1; k <
m_kk; k++) {
3669 molality[j]*molality[k]*
m_CMX_IJ_P[counterIJ2];
3678 for (
size_t k = 1; k <
m_kk; k++) {
3685 if (zeta_P != 0.0) {
3686 sum5 += molality[j]*molality[k]*zeta_P;
3693 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3700 for (
size_t j = 1; j <
m_kk; j++) {
3704 for (
size_t k = 1; k <
m_kk; k++) {
3712 double sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_P[i];
3731 double term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3733 for (
size_t j = 1; j <
m_kk; j++) {
3736 for (
size_t k = 1; k <
m_kk; k++) {
3739 size_t n =
m_kk*j + k;
3741 sum1 += molality[j]*molality[k]*
3746 for (
size_t k = j+1; k <
m_kk; k++) {
3747 if (j == (
m_kk-1)) {
3749 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3750 "logic error 1 in Step 9 of hmw_act");
3755 size_t n =
m_kk*j + k;
3758 for (
size_t m = 1; m <
m_kk; m++) {
3762 sum2 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_P[n];
3771 for (
size_t k = j+1; k <
m_kk; k++) {
3774 throw CanteraError(
"HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3775 "logic error 2 in Step 9 of hmw_act");
3780 size_t n =
m_kk*j + k;
3784 for (
size_t m = 1; m <
m_kk; m++) {
3787 sum3 += molality[j]*molality[k]*molality[m]*
m_Psi_ijk_P[n];
3796 for (
size_t k = 1; k <
m_kk; k++) {
3806 }
else if (k == j) {
3812 for (
size_t m = 1; m <
m_kk; m++) {
3817 if (zeta_P != 0.0) {
3818 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
3825 sum7 += molality[j] * molality[j] * molality[j] *
m_Mu_nnn_P[j];
3828 double sum_m_phi_minus_1 = 2.0 *
3829 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3833 double d_osmotic_coef_dP;
3834 if (molalitysum > 1.0E-150) {
3835 d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3837 d_osmotic_coef_dP = 0.0;
3839 double d_lnwateract_dP = -(
m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
3852 if( m_last_is == is ) {
3859 double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
3860 double aphi = 0.392;
3861 if (is < 1.0E-150) {
3862 for (
int i = 0; i < 17; i++) {
3871 for (
int i=1; i<=4; i++) {
3872 for (
int j=i; j<=4; j++) {
3876 double zprod = (double)ij;
3879 double x = 6.0* zprod * aphi * sqrt(is);
3881 double jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4)));
3883 double t = c3 * c4 * pow(x,c4);
3884 double dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
3885 double jprime = (jfunc/x)*(1.0 + jfunc*dj);
3887 elambda[ij] = zprod*jfunc / (4.0*is);
3888 elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
3895 double* etheta,
double* etheta_prime)
const
3902 "we shouldn't be here");
3904 "called with one species being neutral");
3910 *etheta_prime = 0.0;
3913 double f1 = (double)i / (2.0 * j);
3914 double f2 = (double)j / (2.0 * i);
3929 for (
size_t k = 1; k <
m_kk; k++) {
3936 double eterm = std::exp(-xoverc);
3938 double fptmp = IMS_bfCut_ - IMS_afCut_ /
IMS_cCut_ - IMS_bfCut_*xoverc
3939 + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
3940 double f_prime = 1.0 + eterm*fptmp;
3941 double f = xmolSolvent + IMS_efCut_
3942 + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
3944 double gptmp = IMS_bgCut_ - IMS_agCut_ /
IMS_cCut_ - IMS_bgCut_*xoverc
3945 + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
3946 double g_prime = 1.0 + eterm*gptmp;
3947 double g = xmolSolvent + IMS_egCut_
3948 + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
3950 double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
3951 double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
3952 double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
3954 tmp = log(xx) + lngammak;
3955 for (
size_t k = 1; k <
m_kk; k++) {
3966 vector<double>& moleF =
m_tmpV;
3972 writelog(
"Index Name MoleF MolalityCropped Charge\n");
3973 for (
size_t k = 0; k <
m_kk; k++) {
3974 writelogf(
"%2d %-16s %14.7le %14.7le %5.1f \n",
3978 writelog(
"\n Species Species beta0MX "
3979 "beta1MX beta2MX CphiMX alphaMX thetaij\n");
3980 for (
size_t i = 1; i <
m_kk - 1; i++) {
3981 for (
size_t j = i+1; j <
m_kk; j++) {
3982 size_t n = i *
m_kk + j;
3984 writelogf(
" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
3992 writelog(
"\n Species Species Species psi \n");
3993 for (
size_t i = 1; i <
m_kk; i++) {
3994 for (
size_t j = 1; j <
m_kk; j++) {
3995 for (
size_t k = 1; k <
m_kk; k++) {
3998 writelogf(
" %-16s %-16s %-16s %9.5f \n",
4015 double afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4016 for (
size_t k = 0; k <
m_kk; k++) {
4017 acMolality[k] *= exp(
charge(k) * afac);
4030 double afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4031 for (
size_t k = 0; k <
m_kk; k++) {
4045 double afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
4046 for (
size_t k = 0; k <
m_kk; k++) {
4060 double afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
4061 for (
size_t k = 0; k <
m_kk; k++) {
4075 double afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
4076 for (
size_t k = 0; k <
m_kk; k++) {
4085 double lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
4086 return lnGammaClMs2;
4093 return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
4100 return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
4107 return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
Headers for the HMWSoln ThermoPhase object, which models concentrated electrolyte solutions (see Ther...
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
A map of string keys to values whose type can vary at runtime.
size_t size() const
Returns the number of elements in this map.
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Base class for exceptions thrown by Cantera classes.
void calcIMSCutoffParams_()
Precalculate the IMS Cutoff parameters for typeCutoff = 2.
vector< double > m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT.
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
vector< double > m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
vector< double > m_Beta1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
void applyphScale(double *acMolality) const override
Apply the current phScale to a set of activity Coefficients or activities.
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer's papers.
double m_A_Debye
A_Debye: this expression appears on the top of the ln actCoeff term in the general Debye-Huckel expre...
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...
vector< double > m_h2func_IJ
hfunc2, was called gprime in Pitzer's paper.
vector< double > m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT. Vector index is counterIJ.
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
vector< double > m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
double IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature.
vector< double > m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
void calcMCCutoffParams_()
Calculate molality cut-off parameters.
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer's papers.
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
vector< double > m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P. Vector index is counterIJ.
vector< double > m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT. Vector index is counterIJ.
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
vector< double > m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
HMWSoln(const string &inputFile="", const string &id="")
Construct and initialize an HMWSoln ThermoPhase object directly from an input file.
vector< double > m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
vector< double > m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P. Vector index is counterIJ.
vector< double > m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T. Vector index is counterIJ.
vector< double > m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT. Vector index is counterIJ.
vector< double > m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T. Vector index is counterIJ.
vector< double > m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT. Vector index is counterIJ.
vector< double > m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT. Vector index is counterIJ.
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
double s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale.
vector< double > m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure.
vector< double > m_tmpV
vector of size m_kk, used as a temporary holding area.
void printCoeffs() const
Print out all of the input Pitzer coefficients.
void getActivityConcentrations(double *c) const override
This method returns an array of generalized activity concentrations.
vector< double > m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
vector< double > m_hfunc_IJ
hfunc, was called gprime in Pitzer's paper.
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
void s_update_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
vector< double > m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P. Vector index is counterIJ.
vector< int > CROP_speciesCropped_
This is a boolean-type vector indicating whether a species's activity coefficient is in the cropped r...
vector< double > m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_V.
vector< double > m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T. Vector index is counterIJ.
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
vector< double > m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P. Vector index is counterIJ.
vector< double > m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P. Vector index is counterIJ.
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
vector< double > m_BMX_IJ
Intermediate variable called BMX in Pitzer's paper.
vector< double > m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P.
vector< double > m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT. Vector index is counterIJ.
vector< double > m_BMX_IJ_P
Derivative of BMX_IJ wrt P. Vector index is counterIJ.
double cv_mole() const override
Molar heat capacity at constant volume. Units: J/kmol/K.
vector< double > m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
void getUnscaledMolalityActivityCoefficients(double *acMolality) const override
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
void setA_Debye(double A)
Set the A_Debye parameter.
vector< double > IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model.
vector< double > m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T.
vector< double > m_BMX_IJ_L
Derivative of BMX_IJ wrt T. Vector index is counterIJ.
virtual double relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
vector< double > m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T. Vector index is counterIJ.
vector< double > m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P. Vector index is counterIJ.
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients' dependence on the temperature.
double s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
vector< double > m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
vector< double > m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T. Vector index is counterIJ.
vector< double > m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients.
vector< double > m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer's paper.
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
vector< double > m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer's paper.
vector< double > m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P. Vector index is counterIJ.
vector< double > m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer's paper.
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
vector< int > m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species.
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
virtual double relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_L.
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
PDSS * m_waterSS
Water standard state calculator.
void calcDensity() override
Calculate the density of the mixture using the partial molar volumes and mole fractions as input.
void calc_lambdas(double is) const
Calculate the lambda interactions.
double elambda1[17]
This is elambda1, MEC.
vector< double > m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT. Vector index is counterIJ.
vector< double > m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
double s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
vector< double > m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T.
vector< double > m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
vector< double > m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT. Vector index is counterIJ.
vector< double > m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer's paper.
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients.
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
vector< double > m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P. Vector index is counterIJ.
void getActivities(double *ac) const override
Get the array of non-dimensional activities at the current solution temperature, pressure,...
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...
vector< double > m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T. Vector index is counterIJ.
vector< double > m_Phi_IJ
Intermediate variable called Phi in Pitzer's paper.
vector< double > m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT. Vector index is counterIJ.
vector< double > m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT.
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
double IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process;
double IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
vector< double > m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT. Vector index is counterIJ.
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
vector< double > m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T. Vector index is counterIJ.
void getPartialMolarCp(double *cpbar) const override
Return an array of partial molar heat capacities for the species in the mixture.
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
vector< double > m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
void initLengths()
Initialize all of the species-dependent lengths in the object.
vector< double > m_g2func_IJ
This is the value of g2(x2) in Pitzer's papers. Vector index is counterIJ.
double gibbs_mole() const override
Molar Gibbs function. Units: J/kmol.
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
vector< double > m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
vector< double > m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P. Vector index is counterIJ.
vector< double > m_CMX_IJ
Intermediate variable called CMX in Pitzer's paper.
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
double MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process;
double satPressure(double T) override
Get the saturation pressure for a given temperature.
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
vector< double > m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT. Vector index is counterIJ.
double elambda[17]
This is elambda, MEC.
vector< double > m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P.
vector< double > m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
vector< double > m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
void calcMolalitiesCropped() const
Calculate the cropped molalities.
vector< double > m_CphiMX_ij_L
Derivative of Cphi_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.
void s_updateScaling_pHScaling_dT2() const
Apply the current phScale to a set of 2nd derivatives of the activity Coefficients wrt temperature.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer's papers.
vector< double > m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P. Vector index is counterIJ.
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
vector< double > m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T. Vector index is counterIJ.
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
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 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...
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
vector< double > m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T. Vector index is counterIJ.
vector< double > m_gamma_tmp
Intermediate storage of the activity coefficient itself.
double s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale.
double m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
size_t m_indexCLM
Index of the phScale species.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
double m_xmolSolventMIN
In any molality implementation, it makes sense to have a minimum solvent mole fraction requirement,...
vector< double > m_molalities
Current value of the molalities of the species in the phase.
int m_pHScalingType
Scaling to be used for output of single-ion species activity coefficients.
void setMoleFSolventMin(double xmolSolventMIN)
Sets the minimum mole fraction in the molality formulation.
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
double m_weightSolvent
Molecular weight of the Solvent.
Class for the liquid water pressure dependent standard state.
virtual double satPressure(double T)
saturation pressure
virtual void setState_TP(double temp, double pres)
Set the internal temperature and pressure.
void assignDensity(const double density_)
Set the internally stored constant density (kg/m^3) of the phase.
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
size_t m_kk
Number of species in the phase.
double temperature() const
Temperature (K).
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
string speciesName(size_t k) const
Name of the species with index k.
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
double moleFraction(size_t k) const
Return the mole fraction of a single species.
int stateMFNumber() const
Return the State Mole Fraction Number.
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
virtual double molarVolume() const
Molar volume (m^3/kmol).
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
virtual double thermalExpansionCoeff() const
Return the volumetric thermal expansion coefficient. Units: 1/K.
virtual void getParameters(AnyMap &phaseNode) const
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
double RT() const
Return the Gas Constant multiplied by the current temperature.
virtual double isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
void initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
AnyMap m_input
Data supplied via setParameters.
double pressure() const override
Returns the current pressure of the phase.
void getEntropy_R(double *sr) const override
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
void getStandardChemPotentials(double *mu) const override
Get the array of chemical potentials at unit activity for the species at their standard states at the...
void getCp_R(double *cpr) const override
Get the nondimensional Heat Capacities at constant pressure for the species standard states at the cu...
void getEnthalpy_RT(double *hrt) const override
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
void getStandardVolumes(double *vol) const override
Get the molar volumes of the species standard states at the current T and P of the solution.
virtual void updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
CachedScalar getScalar(int id)
Get a reference to a CachedValue object representing a scalar (double) with the given id.
int getId()
Get a unique id for a cached value.
Header file for a common definitions used in electrolytes thermodynamics.
bool caseInsensitiveEquals(const string &input, const string &test)
Case insensitive equality predicate.
#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.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
int sign(T x)
Sign of a number.
const double GasConstant
Universal Gas Constant [J/kmol/K].
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
const int PHSCALE_PITZER
Scale to be used for the output of single-ion activity coefficients is that used by Pitzer.
const double SmallNumber
smallest number to compare to zero.
const int PHSCALE_NBS
Scale to be used for evaluation of single-ion activity coefficients is that used by the NBS standard ...
Contains declarations for string manipulation functions within Cantera.
A cached property value and the state at which it was evaluated.
T value
The value of the cached property.
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.