Cantera  3.1.0a1
Loading...
Searching...
No Matches
HMWSoln.cpp
Go to the documentation of this file.
1/**
2 * @file HMWSoln.cpp
3 * Definitions for the HMWSoln ThermoPhase object, which
4 * models concentrated electrolyte solutions
5 * (see @ref thermoprops and @link Cantera::HMWSoln HMWSoln @endlink) .
6 *
7 * Class HMWSoln represents a concentrated liquid electrolyte phase which obeys
8 * the Pitzer formulation for nonideality using molality-based standard states.
9 *
10 * This version of the code was modified to have the binary Beta2 Pitzer
11 * parameter consistent with the temperature expansions used for Beta0,
12 * Beta1, and Cphi.(CFJC, SNL)
13 */
14
15// This file is part of Cantera. See License.txt in the top-level directory or
16// at https://cantera.org/license.txt for license and copyright information.
17
23
24namespace Cantera
25{
26
27namespace {
28double A_Debye_default = 1.172576; // units = sqrt(kg/gmol)
29double maxIionicStrength_default = 100.0;
30double crop_ln_gamma_o_min_default = -6.0;
31double crop_ln_gamma_o_max_default = 3.0;
32double crop_ln_gamma_k_min_default = -5.0;
33double crop_ln_gamma_k_max_default = 15.0;
34}
35
36HMWSoln::~HMWSoln()
37{
38 // Defined in .cpp to limit dependence on WaterProps.h
39}
40
41HMWSoln::HMWSoln(const string& inputFile, const string& id_) :
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),
47 m_last_is(-1.0)
48{
49 initThermoFile(inputFile, id_);
50}
51
52// -------- Molar Thermodynamic Properties of the Solution ---------------
53
55{
57 return mean_X(m_tmpV);
58}
59
61{
63 double hbar = mean_X(m_tmpV);
65 for (size_t k = 0; k < m_kk; k++) {
66 m_gamma_tmp[k] *= RT();
67 }
68 double h0bar = mean_X(m_gamma_tmp);
69 return hbar - h0bar;
70}
71
73{
74 double L = relative_enthalpy();
76 double xanion = 0.0;
77 size_t kcation = npos;
78 double xcation = 0.0;
79 size_t kanion = npos;
80 for (size_t k = 0; k < m_kk; k++) {
81 if (charge(k) > 0.0) {
82 if (m_tmpV[k] > xanion) {
83 xanion = m_tmpV[k];
84 kanion = k;
85 }
86 } else if (charge(k) < 0.0) {
87 if (m_tmpV[k] > xcation) {
88 xcation = m_tmpV[k];
89 kcation = k;
90 }
91 }
92 }
93 if (kcation == npos || kanion == npos) {
94 return L;
95 }
96 double xuse = xcation;
97 double factor = 1;
98 if (xanion < xcation) {
99 xuse = xanion;
100 if (charge(kcation) != 1.0) {
101 factor = charge(kcation);
102 }
103 } else {
104 if (charge(kanion) != 1.0) {
105 factor = charge(kanion);
106 }
107 }
108 xuse = xuse / factor;
109 return L / xuse;
110}
111
113{
115 return mean_X(m_tmpV);
116}
117
119{
121 return mean_X(m_tmpV);
122}
123
124double HMWSoln::cp_mole() const
125{
127 return mean_X(m_tmpV);
128}
129
130double HMWSoln::cv_mole() const
131{
132 double kappa_t = isothermalCompressibility();
133 double beta = thermalExpansionCoeff();
134 double cp = cp_mole();
135 double tt = temperature();
136 double molarV = molarVolume();
137 return cp - beta * beta * tt * molarV / kappa_t;
138}
139
140// ------- Mechanical Equation of State Properties ------------------------
141
143{
144 static const int cacheId = m_cache.getId();
145 CachedScalar cached = m_cache.getScalar(cacheId);
146 if(cached.validate(temperature(), pressure(), stateMFNumber())) {
147 return;
148 }
149
150 // Calculate all of the other standard volumes. Note these are constant for
151 // now
153 double dd = meanMolecularWeight() / mean_X(m_tmpV);
155}
156
157// ------- Activities and Activity Concentrations
158
160{
161 double cs_solvent = standardConcentration();
162 getActivities(c);
163 c[0] *= cs_solvent;
164 if (m_kk > 1) {
165 double cs_solute = standardConcentration(1);
166 for (size_t k = 1; k < m_kk; k++) {
167 c[k] *= cs_solute;
168 }
169 }
170}
171
172double HMWSoln::standardConcentration(size_t k) const
173{
175 double mvSolvent = m_tmpV[0];
176 if (k > 0) {
177 return m_Mnaught / mvSolvent;
178 }
179 return 1.0 / mvSolvent;
180}
181
182void HMWSoln::getActivities(double* ac) const
183{
185
186 // Update the molality array, m_molalities(). This requires an update due to
187 // mole fractions
189
190 // Now calculate the array of activities.
191 for (size_t k = 1; k < m_kk; k++) {
192 ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
193 }
194 double xmolSolvent = moleFraction(0);
195 ac[0] = exp(m_lnActCoeffMolal_Scaled[0]) * xmolSolvent;
196}
197
199{
201 A_Debye_TP(-1.0, -1.0);
203 std::copy(m_lnActCoeffMolal_Unscaled.begin(), m_lnActCoeffMolal_Unscaled.end(), acMolality);
204 for (size_t k = 0; k < m_kk; k++) {
205 acMolality[k] = exp(acMolality[k]);
206 }
207}
208
209// ------ Partial Molar Properties of the Solution -----------------
210
211void HMWSoln::getChemPotentials(double* mu) const
212{
213 double xx;
214
215 // First get the standard chemical potentials in molar form. This requires
216 // updates of standard state as a function of T and P
218
219 // Update the activity coefficients. This also updates the internal molality
220 // array.
222 double xmolSolvent = moleFraction(0);
223 for (size_t k = 1; k < m_kk; k++) {
224 xx = std::max(m_molalities[k], SmallNumber);
225 mu[k] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
226 }
227 xx = std::max(xmolSolvent, SmallNumber);
228 mu[0] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[0]);
229}
230
232{
233 // Get the nondimensional standard state enthalpies
234 getEnthalpy_RT(hbar);
235
236 // dimensionalize it.
237 for (size_t k = 0; k < m_kk; k++) {
238 hbar[k] *= RT();
239 }
240
241 // Update the activity coefficients, This also update the internally stored
242 // molalities.
245 for (size_t k = 0; k < m_kk; k++) {
246 hbar[k] -= RT() * temperature() * m_dlnActCoeffMolaldT_Scaled[k];
247 }
248}
249
250void HMWSoln::getPartialMolarEntropies(double* sbar) const
251{
252 // Get the standard state entropies at the temperature and pressure of the
253 // solution.
254 getEntropy_R(sbar);
255
256 // Dimensionalize the entropies
257 for (size_t k = 0; k < m_kk; k++) {
258 sbar[k] *= GasConstant;
259 }
260
261 // Update the activity coefficients, This also update the internally stored
262 // molalities.
264
265 // First we will add in the obvious dependence on the T term out front of
266 // the log activity term
267 double mm;
268 for (size_t k = 1; k < m_kk; k++) {
269 mm = std::max(SmallNumber, m_molalities[k]);
270 sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
271 }
272 double xmolSolvent = moleFraction(0);
273 mm = std::max(SmallNumber, xmolSolvent);
274 sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[0]);
275
276 // Check to see whether activity coefficients are temperature dependent. If
277 // they are, then calculate the their temperature derivatives and add them
278 // into the result.
280 for (size_t k = 0; k < m_kk; k++) {
281 sbar[k] -= RT() * m_dlnActCoeffMolaldT_Scaled[k];
282 }
283}
284
285void HMWSoln::getPartialMolarVolumes(double* vbar) const
286{
287 // Get the standard state values in m^3 kmol-1
288 getStandardVolumes(vbar);
289
290 // Update the derivatives wrt the activity coefficients.
293 for (size_t k = 0; k < m_kk; k++) {
294 vbar[k] += RT() * m_dlnActCoeffMolaldP_Scaled[k];
295 }
296}
297
298void HMWSoln::getPartialMolarCp(double* cpbar) const
299{
300 getCp_R(cpbar);
301 for (size_t k = 0; k < m_kk; k++) {
302 cpbar[k] *= GasConstant;
303 }
304
305 // Update the activity coefficients, This also update the internally stored
306 // molalities.
310 for (size_t k = 0; k < m_kk; k++) {
311 cpbar[k] -= (2.0 * RT() * m_dlnActCoeffMolaldT_Scaled[k] +
313 }
314}
315
316// -------------- Utilities -------------------------------
317
318double HMWSoln::satPressure(double t) {
319 double p_old = pressure();
320 double t_old = temperature();
321 double pres = m_waterSS->satPressure(t);
322
323 // Set the underlying object back to its original state.
324 m_waterSS->setState_TP(t_old, p_old);
325 return pres;
326}
327
328static void check_nParams(const string& method, size_t nParams, size_t m_formPitzerTemp)
329{
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);
336 }
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);
340 }
341}
342
343void 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)
346{
347 size_t k1 = speciesIndex(sp1);
348 size_t k2 = speciesIndex(sp2);
349 if (k1 == npos) {
350 throw CanteraError("HMWSoln::setBinarySalt", "Species '{}' not found", sp1);
351 } else if (k2 == npos) {
352 throw CanteraError("HMWSoln::setBinarySalt", "Species '{}' not found", sp2);
353 }
354 if (charge(k1) < 0 && charge(k2) > 0) {
355 std::swap(k1, k2);
356 } else if (charge(k1) * charge(k2) >= 0) {
357 throw CanteraError("HMWSoln::setBinarySalt", "Species '{}' and '{}' "
358 "do not have opposite charges ({}, {})", sp1, sp2,
359 charge(k1), charge(k2));
360 }
361 check_nParams("HMWSoln::setBinarySalt", nParams, m_formPitzerTemp);
362
363 size_t c = m_CounterIJ[k1 * m_kk + k2];
364 m_Beta0MX_ij[c] = beta0[0];
365 m_Beta1MX_ij[c] = beta1[0];
366 m_Beta2MX_ij[c] = beta2[0];
367 m_CphiMX_ij[c] = Cphi[0];
368 for (size_t n = 0; n < nParams; n++) {
369 m_Beta0MX_ij_coeff(n, c) = beta0[n];
370 m_Beta1MX_ij_coeff(n, c) = beta1[n];
371 m_Beta2MX_ij_coeff(n, c) = beta2[n];
372 m_CphiMX_ij_coeff(n, c) = Cphi[n];
373 }
374 m_Alpha1MX_ij[c] = alpha1;
375 m_Alpha2MX_ij[c] = alpha2;
376}
377
378void HMWSoln::setTheta(const string& sp1, const string& sp2,
379 size_t nParams, double* theta)
380{
381 size_t k1 = speciesIndex(sp1);
382 size_t k2 = speciesIndex(sp2);
383 if (k1 == npos) {
384 throw CanteraError("HMWSoln::setTheta", "Species '{}' not found", sp1);
385 } else if (k2 == npos) {
386 throw CanteraError("HMWSoln::setTheta", "Species '{}' not found", sp2);
387 }
388 if (charge(k1) * charge(k2) <= 0) {
389 throw CanteraError("HMWSoln::setTheta", "Species '{}' and '{}' "
390 "should both have the same (non-zero) charge ({}, {})", sp1, sp2,
391 charge(k1), charge(k2));
392 }
393 check_nParams("HMWSoln::setTheta", nParams, m_formPitzerTemp);
394 size_t c = m_CounterIJ[k1 * m_kk + k2];
395 m_Theta_ij[c] = theta[0];
396 for (size_t n = 0; n < nParams; n++) {
397 m_Theta_ij_coeff(n, c) = theta[n];
398 }
399}
400
401void HMWSoln::setPsi(const string& sp1, const string& sp2,
402 const string& sp3, size_t nParams, double* psi)
403{
404 size_t k1 = speciesIndex(sp1);
405 size_t k2 = speciesIndex(sp2);
406 size_t k3 = speciesIndex(sp3);
407 if (k1 == npos) {
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);
413 }
414
415 if (!charge(k1) || !charge(k2) || !charge(k3) ||
416 std::abs(sign(charge(k1) + sign(charge(k2)) + sign(charge(k3)))) != 1) {
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 {} ({}).",
420 sp1, charge(k1), sp2, charge(k2), sp3, charge(k3));
421 }
422 check_nParams("HMWSoln::setPsi", nParams, m_formPitzerTemp);
423 auto cc = {k1*m_kk*m_kk + k2*m_kk + k3,
424 k1*m_kk*m_kk + k3*m_kk + k2,
425 k2*m_kk*m_kk + k1*m_kk + k3,
426 k2*m_kk*m_kk + k3*m_kk + k1,
427 k3*m_kk*m_kk + k2*m_kk + k1,
428 k3*m_kk*m_kk + k1*m_kk + k2};
429 for (auto c : cc) {
430 for (size_t n = 0; n < nParams; n++) {
431 m_Psi_ijk_coeff(n, c) = psi[n];
432 }
433 m_Psi_ijk[c] = psi[0];
434 }
435}
436
437void HMWSoln::setLambda(const string& sp1, const string& sp2,
438 size_t nParams, double* lambda)
439{
440 size_t k1 = speciesIndex(sp1);
441 size_t k2 = speciesIndex(sp2);
442 if (k1 == npos) {
443 throw CanteraError("HMWSoln::setLambda", "Species '{}' not found", sp1);
444 } else if (k2 == npos) {
445 throw CanteraError("HMWSoln::setLambda", "Species '{}' not found", sp2);
446 }
447
448 if (charge(k1) != 0 && charge(k2) != 0) {
449 throw CanteraError("HMWSoln::setLambda", "Expected at least one neutral"
450 " species, but given species (charges) were: {} ({}) and {} ({}).",
451 sp1, charge(k1), sp2, charge(k2));
452 }
453 if (charge(k1) != 0) {
454 std::swap(k1, k2);
455 }
456 check_nParams("HMWSoln::setLambda", nParams, m_formPitzerTemp);
457 size_t c = k1*m_kk + k2;
458 for (size_t n = 0; n < nParams; n++) {
459 m_Lambda_nj_coeff(n, c) = lambda[n];
460 }
461 m_Lambda_nj(k1, k2) = lambda[0];
462}
463
464void HMWSoln::setMunnn(const string& sp, size_t nParams, double* munnn)
465{
466 size_t k = speciesIndex(sp);
467 if (k == npos) {
468 throw CanteraError("HMWSoln::setMunnn", "Species '{}' not found", sp);
469 }
470
471 if (charge(k) != 0) {
472 throw CanteraError("HMWSoln::setMunnn", "Expected a neutral species,"
473 " got {} ({}).", sp, charge(k));
474 }
475 check_nParams("HMWSoln::setMunnn", nParams, m_formPitzerTemp);
476 for (size_t n = 0; n < nParams; n++) {
477 m_Mu_nnn_coeff(n, k) = munnn[n];
478 }
479 m_Mu_nnn[k] = munnn[0];
480}
481
482void HMWSoln::setZeta(const string& sp1, const string& sp2,
483 const string& sp3, size_t nParams, double* psi)
484{
485 size_t k1 = speciesIndex(sp1);
486 size_t k2 = speciesIndex(sp2);
487 size_t k3 = speciesIndex(sp3);
488 if (k1 == npos) {
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);
494 }
495
496 if (charge(k1)*charge(k2)*charge(k3) != 0 ||
497 sign(charge(k1)) + sign(charge(k2)) + sign(charge(k3)) != 0) {
498 throw CanteraError("HMWSoln::setZeta", "Requires one neutral species, "
499 "one cation, and one anion, but given species (charges) were: "
500 "{} ({}), {} ({}), and {} ({}).",
501 sp1, charge(k1), sp2, charge(k2), sp3, charge(k3));
502 }
503
504 //! Make k1 the neutral species
505 if (charge(k2) == 0) {
506 std::swap(k1, k2);
507 } else if (charge(k3) == 0) {
508 std::swap(k1, k3);
509 }
510
511 // Make k2 the cation
512 if (charge(k3) > 0) {
513 std::swap(k2, k3);
514 }
515
516 check_nParams("HMWSoln::setZeta", nParams, m_formPitzerTemp);
517 // In contrast to setPsi, there are no duplicate entries
518 size_t c = k1 * m_kk *m_kk + k2 * m_kk + k3;
519 for (size_t n = 0; n < nParams; n++) {
520 m_Psi_ijk_coeff(n, c) = psi[n];
521 }
522 m_Psi_ijk[c] = psi[0];
523}
524
525void HMWSoln::setPitzerTempModel(const string& model)
526{
527 if (caseInsensitiveEquals(model, "constant") || caseInsensitiveEquals(model, "default")) {
528 m_formPitzerTemp = PITZER_TEMP_CONSTANT;
529 } else if (caseInsensitiveEquals(model, "linear")) {
530 m_formPitzerTemp = PITZER_TEMP_LINEAR;
531 } else if (caseInsensitiveEquals(model, "complex") || caseInsensitiveEquals(model, "complex1")) {
532 m_formPitzerTemp = PITZER_TEMP_COMPLEX1;
533 } else {
534 throw CanteraError("HMWSoln::setPitzerTempModel",
535 "Unknown Pitzer ActivityCoeff Temp model: {}", model);
536 }
537}
538
540{
541 if (A < 0) {
542 m_form_A_Debye = A_DEBYE_WATER;
543 } else {
544 m_form_A_Debye = A_DEBYE_CONST;
545 m_A_Debye = A;
546 }
547}
548
549void HMWSoln::setCroppingCoefficients(double ln_gamma_k_min,
550 double ln_gamma_k_max, double ln_gamma_o_min, double ln_gamma_o_max)
551{
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;
556}
557
558vector<double> getSizedVector(const AnyMap& item, const string& key, size_t nCoeffs)
559{
560 vector<double> v;
561 if (item[key].is<double>()) {
562 // Allow a single value to be given directly, rather than as a list of
563 // one item
564 v.push_back(item[key].asDouble());
565 } else {
566 v = item[key].asVector<double>(1, nCoeffs);
567 }
568 if (v.size() == 1 && nCoeffs == 5) {
569 // Adapt constant-temperature data to be compatible with the "complex"
570 // temperature model
571 v.resize(5, 0.0);
572 }
573 return v;
574}
575
577{
579 if (m_input.hasKey("activity-data")) {
580 auto& actData = m_input["activity-data"].as<AnyMap>();
581 setPitzerTempModel(actData["temperature-model"].asString());
582 initLengths();
583 size_t nCoeffs = 1;
584 if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
585 nCoeffs = 2;
586 } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
587 nCoeffs = 5;
588 }
589 if (actData.hasKey("A_Debye")) {
590 if (actData["A_Debye"] == "variable") {
591 setA_Debye(-1);
592 } else {
593 setA_Debye(actData.convert("A_Debye", "kg^0.5/gmol^0.5"));
594 }
595 }
596 if (actData.hasKey("max-ionic-strength")) {
597 setMaxIonicStrength(actData["max-ionic-strength"].asDouble());
598 }
599 if (actData.hasKey("interactions")) {
600 for (auto& item : actData["interactions"].asVector<AnyMap>()) {
601 auto& species = item["species"].asVector<string>(1, 3);
602 size_t nsp = species.size();
603 double q0 = charge(speciesIndex(species[0]));
604 double q1 = (nsp > 1) ? charge(speciesIndex(species[1])) : 0;
605 double q2 = (nsp == 3) ? charge(speciesIndex(species[2])) : 0;
606 if (nsp == 2 && q0 * q1 < 0) {
607 // Two species with opposite charges - binary salt
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()) {
614 throw InputFileError("HMWSoln::initThermo", item,
615 "Inconsistent binary salt array sizes ({}, {}, {}, {})",
616 beta0.size(), beta1.size(), beta2.size(), Cphi.size());
617 }
618 double alpha1 = item["alpha1"].asDouble();
619 double alpha2 = item.getDouble("alpha2", 0.0);
620 setBinarySalt(species[0], species[1], beta0.size(),
621 beta0.data(), beta1.data(), beta2.data(), Cphi.data(),
622 alpha1, alpha2);
623 } else if (nsp == 2 && q0 * q1 > 0) {
624 // Two species with like charges - "theta" interaction
625 vector<double> theta = getSizedVector(item, "theta", nCoeffs);
626 setTheta(species[0], species[1], theta.size(), theta.data());
627 } else if (nsp == 2 && q0 * q1 == 0) {
628 // Two species, including at least one neutral
629 vector<double> lambda = getSizedVector(item, "lambda", nCoeffs);
630 setLambda(species[0], species[1], lambda.size(), lambda.data());
631 } else if (nsp == 3 && q0 * q1 * q2 != 0) {
632 // Three charged species - "psi" interaction
633 vector<double> psi = getSizedVector(item, "psi", nCoeffs);
634 setPsi(species[0], species[1], species[2],
635 psi.size(), psi.data());
636 } else if (nsp == 3 && q0 * q1 * q2 == 0) {
637 // Three species, including one neutral
638 vector<double> zeta = getSizedVector(item, "zeta", nCoeffs);
639 setZeta(species[0], species[1], species[2],
640 zeta.size(), zeta.data());
641 } else if (nsp == 1) {
642 // single species (should be neutral)
643 vector<double> mu = getSizedVector(item, "mu", nCoeffs);
644 setMunnn(species[0], mu.size(), mu.data());
645 }
646 }
647 }
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));
655 }
656 } else {
657 initLengths();
658 }
659
660 for (int i = 0; i < 17; i++) {
661 elambda[i] = 0.0;
662 elambda1[i] = 0.0;
663 }
664
665 // Store a local pointer to the water standard state model.
666 m_waterSS = providePDSS(0);
667
668 // Initialize the water property calculator. It will share the internal eos
669 // water calculator.
670 m_waterProps = make_unique<WaterProps>(dynamic_cast<PDSS_Water*>(m_waterSS));
671
672 // Lastly calculate the charge balance and then add stuff until the charges
673 // compensate
674 vector<double> mf(m_kk, 0.0);
675 getMoleFractions(mf.data());
676 bool notDone = true;
677
678 while (notDone) {
679 double sum = 0.0;
680 size_t kMaxC = npos;
681 double MaxC = 0.0;
682 for (size_t k = 0; k < m_kk; k++) {
683 sum += mf[k] * charge(k);
684 if (fabs(mf[k] * charge(k)) > MaxC) {
685 kMaxC = k;
686 }
687 }
688 size_t kHp = speciesIndex("H+");
689 size_t kOHm = speciesIndex("OH-");
690
691 if (fabs(sum) > 1.0E-30) {
692 if (kHp != npos) {
693 if (mf[kHp] > sum * 1.1) {
694 mf[kHp] -= sum;
695 mf[0] += sum;
696 notDone = false;
697 } else {
698 if (sum > 0.0) {
699 mf[kHp] *= 0.5;
700 mf[0] += mf[kHp];
701 sum -= mf[kHp];
702 }
703 }
704 }
705 if (notDone) {
706 if (kOHm != npos) {
707 if (mf[kOHm] > -sum * 1.1) {
708 mf[kOHm] += sum;
709 mf[0] -= sum;
710 notDone = false;
711 } else {
712 if (sum < 0.0) {
713 mf[kOHm] *= 0.5;
714 mf[0] += mf[kOHm];
715 sum += mf[kOHm];
716 }
717 }
718 }
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);
723 } else {
724 mf[kMaxC] *= 0.5;
725 mf[0] += mf[kMaxC];
726 notDone = true;
727 }
728 }
729 }
730 setMoleFractions(mf.data());
731 } else {
732 notDone = false;
733 }
734 }
735
738 setMoleFSolventMin(1.0E-5);
739}
740
741void assignTrimmed(AnyMap& interaction, const string& key, vector<double>& values) {
742 while (values.size() > 1 && values.back() == 0) {
743 values.pop_back();
744 }
745 if (values.size() == 1) {
746 interaction[key] = values[0];
747 } else {
748 interaction[key] = values;
749 }
750}
751
752void HMWSoln::getParameters(AnyMap& phaseNode) const
753{
755 AnyMap activityNode;
756 size_t nParams = 1;
757 if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
758 activityNode["temperature-model"] = "linear";
759 nParams = 2;
760 } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
761 activityNode["temperature-model"] = "complex";
762 nParams = 5;
763 }
764
765 if (m_form_A_Debye == A_DEBYE_WATER) {
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");
769 }
770 if (m_maxIionicStrength != maxIionicStrength_default) {
771 activityNode["max-ionic-strength"] = m_maxIionicStrength;
772 }
773
774 vector<AnyMap> interactions;
775
776 // Binary 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;
780 // lambda: neutral-charged / neutral-neutral interactions
781 bool lambda_found = false;
782 for (size_t n = 0; n < nParams; n++) {
783 if (m_Lambda_nj_coeff(n, c)) {
784 lambda_found = true;
785 break;
786 }
787 }
788 if (lambda_found) {
789 AnyMap interaction;
790 interaction["species"] = vector<string>{
791 speciesName(i), speciesName(j)};
792 vector<double> lambda(nParams);
793 for (size_t n = 0; n < nParams; n++) {
794 lambda[n] = m_Lambda_nj_coeff(n, c);
795 }
796 assignTrimmed(interaction, "lambda", lambda);
797 interactions.push_back(std::move(interaction));
798 continue;
799 }
800
801 c = static_cast<size_t>(m_CounterIJ[m_kk * i + j]);
802 if (c == 0 || i > j) {
803 continue;
804 }
805
806 // beta: opposite charged interactions
807 bool salt_found = false;
808 for (size_t n = 0; n < nParams; n++) {
809 if (m_Beta0MX_ij_coeff(n, c) || m_Beta1MX_ij_coeff(n, c) ||
811 {
812 salt_found = true;
813 break;
814 }
815 }
816 if (salt_found) {
817 AnyMap interaction;
818 interaction["species"] = vector<string>{
819 speciesName(i), speciesName(j)};
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++) {
823 beta0[n] = m_Beta0MX_ij_coeff(n, c);
824 beta1[n] = m_Beta1MX_ij_coeff(n, c);
825 beta2[n] = m_Beta2MX_ij_coeff(n, c);
826 Cphi[n] = m_CphiMX_ij_coeff(n, c);
827 if (beta0[n] || beta1[n] || beta2[n] || Cphi[n]) {
828 last_nonzero = n;
829 }
830 }
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];
836 } else {
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;
845 }
846 interaction["alpha1"] = m_Alpha1MX_ij[c];
847 if (m_Alpha2MX_ij[c]) {
848 interaction["alpha2"] = m_Alpha2MX_ij[c];
849 }
850 interactions.push_back(std::move(interaction));
851 continue;
852 }
853
854 // theta: like-charge interactions
855 bool theta_found = false;
856 for (size_t n = 0; n < nParams; n++) {
857 if (m_Theta_ij_coeff(n, c)) {
858 theta_found = true;
859 break;
860 }
861 }
862 if (theta_found) {
863 AnyMap interaction;
864 interaction["species"] = vector<string>{
865 speciesName(i), speciesName(j)};
866 vector<double> theta(nParams);
867 for (size_t n = 0; n < nParams; n++) {
868 theta[n] = m_Theta_ij_coeff(n, c);
869 }
870 assignTrimmed(interaction, "theta", theta);
871 interactions.push_back(std::move(interaction));
872 continue;
873 }
874 }
875 }
876
877 // psi: ternary charged species interactions
878 // Need to check species charges because both psi and zeta parameters
879 // are stored in m_Psi_ijk_coeff
880 for (size_t i = 1; i < m_kk; i++) {
881 if (charge(i) == 0) {
882 continue;
883 }
884 for (size_t j = i + 1; j < m_kk; j++) {
885 if (charge(j) == 0) {
886 continue;
887 }
888 for (size_t k = j + 1; k < m_kk; k++) {
889 if (charge(k) == 0) {
890 continue;
891 }
892 size_t c = i*m_kk*m_kk + j*m_kk + k;
893 for (size_t n = 0; n < nParams; n++) {
894 if (m_Psi_ijk_coeff(n, c) != 0) {
895 AnyMap interaction;
896 interaction["species"] = vector<string>{
898 vector<double> psi(nParams);
899 for (size_t m = 0; m < nParams; m++) {
900 psi[m] = m_Psi_ijk_coeff(m, c);
901 }
902 assignTrimmed(interaction, "psi", psi);
903 interactions.push_back(std::move(interaction));
904 break;
905 }
906 }
907 }
908 }
909 }
910
911 // zeta: neutral-cation-anion interactions
912 for (size_t i = 1; i < m_kk; i++) {
913 if (charge(i) != 0) {
914 continue; // first species must be neutral
915 }
916 for (size_t j = 1; j < m_kk; j++) {
917 if (charge(j) <= 0) {
918 continue; // second species must be cation
919 }
920 for (size_t k = 1; k < m_kk; k++) {
921 if (charge(k) >= 0) {
922 continue; // third species must be anion
923 }
924 size_t c = i*m_kk*m_kk + j*m_kk + k;
925 for (size_t n = 0; n < nParams; n++) {
926 if (m_Psi_ijk_coeff(n, c) != 0) {
927 AnyMap interaction;
928 interaction["species"] = vector<string>{
930 vector<double> zeta(nParams);
931 for (size_t m = 0; m < nParams; m++) {
932 zeta[m] = m_Psi_ijk_coeff(m, c);
933 }
934 assignTrimmed(interaction, "zeta", zeta);
935 interactions.push_back(std::move(interaction));
936 break;
937 }
938 }
939 }
940 }
941 }
942
943 // mu: neutral self-interaction
944 for (size_t i = 1; i < m_kk; i++) {
945 for (size_t n = 0; n < nParams; n++) {
946 if (m_Mu_nnn_coeff(n, i) != 0) {
947 AnyMap interaction;
948 interaction["species"] = vector<string>{speciesName(i)};
949 vector<double> mu(nParams);
950 for (size_t m = 0; m < nParams; m++) {
951 mu[m] = m_Mu_nnn_coeff(m, i);
952 }
953 assignTrimmed(interaction, "mu", mu);
954 interactions.push_back(std::move(interaction));
955 break;
956 }
957 }
958 }
959
960 activityNode["interactions"] = std::move(interactions);
961
962 AnyMap croppingCoeffs;
963 if (CROP_ln_gamma_k_min != crop_ln_gamma_k_min_default) {
964 croppingCoeffs["ln_gamma_k_min"] = CROP_ln_gamma_k_min;
965 }
966 if (CROP_ln_gamma_k_max != crop_ln_gamma_k_max_default) {
967 croppingCoeffs["ln_gamma_k_max"] = CROP_ln_gamma_k_max;
968 }
969 if (CROP_ln_gamma_o_min != crop_ln_gamma_o_min_default) {
970 croppingCoeffs["ln_gamma_o_min"] = CROP_ln_gamma_o_min;
971 }
972 if (CROP_ln_gamma_o_max != crop_ln_gamma_o_max_default) {
973 croppingCoeffs["ln_gamma_o_max"] = CROP_ln_gamma_o_max;
974 }
975 if (croppingCoeffs.size()) {
976 activityNode["cropping-coefficients"] = std::move(croppingCoeffs);
977 }
978
979 phaseNode["activity-data"] = std::move(activityNode);
980}
981
982double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
983{
984 double T = temperature();
985 double A;
986 if (tempArg != -1.0) {
987 T = tempArg;
988 }
989 double P = pressure();
990 if (presArg != -1.0) {
991 P = presArg;
992 }
993
994 static const int cacheId = m_cache.getId();
995 CachedScalar cached = m_cache.getScalar(cacheId);
996 if(cached.validate(T, P)) {
997 return m_A_Debye;
998 }
999
1000 switch (m_form_A_Debye) {
1001 case A_DEBYE_CONST:
1002 A = m_A_Debye;
1003 break;
1004 case A_DEBYE_WATER:
1005 A = m_waterProps->ADebye(T, P, 0);
1006 m_A_Debye = A;
1007 break;
1008 default:
1009 throw CanteraError("HMWSoln::A_Debye_TP", "shouldn't be here");
1010 }
1011 return A;
1012}
1013
1014double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const
1015{
1016 double T = temperature();
1017 if (tempArg != -1.0) {
1018 T = tempArg;
1019 }
1020 double P = pressure();
1021 if (presArg != -1.0) {
1022 P = presArg;
1023 }
1024 double dAdT;
1025 switch (m_form_A_Debye) {
1026 case A_DEBYE_CONST:
1027 dAdT = 0.0;
1028 break;
1029 case A_DEBYE_WATER:
1030 dAdT = m_waterProps->ADebye(T, P, 1);
1031 break;
1032 default:
1033 throw CanteraError("HMWSoln::dA_DebyedT_TP", "shouldn't be here");
1034 }
1035 return dAdT;
1036}
1037
1038double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
1039{
1040 double T = temperature();
1041 if (tempArg != -1.0) {
1042 T = tempArg;
1043 }
1044 double P = pressure();
1045 if (presArg != -1.0) {
1046 P = presArg;
1047 }
1048
1049 double dAdP;
1050 static const int cacheId = m_cache.getId();
1051 CachedScalar cached = m_cache.getScalar(cacheId);
1052 switch (m_form_A_Debye) {
1053 case A_DEBYE_CONST:
1054 dAdP = 0.0;
1055 break;
1056 case A_DEBYE_WATER:
1057 if(cached.validate(T, P)) {
1058 dAdP = cached.value;
1059 } else {
1060 dAdP = m_waterProps->ADebye(T, P, 3);
1061 cached.value = dAdP;
1062 }
1063 break;
1064 default:
1065 throw CanteraError("HMWSoln::dA_DebyedP_TP", "shouldn't be here");
1066 }
1067 return dAdP;
1068}
1069
1070double HMWSoln::ADebye_L(double tempArg, double presArg) const
1071{
1072 double dAdT = dA_DebyedT_TP();
1073 double dAphidT = dAdT /3.0;
1074 double T = temperature();
1075 if (tempArg != -1.0) {
1076 T = tempArg;
1077 }
1078 return dAphidT * (4.0 * GasConstant * T * T);
1079}
1080
1081double HMWSoln::ADebye_V(double tempArg, double presArg) const
1082{
1083 double dAdP = dA_DebyedP_TP();
1084 double dAphidP = dAdP /3.0;
1085 double T = temperature();
1086 if (tempArg != -1.0) {
1087 T = tempArg;
1088 }
1089 return - dAphidP * (4.0 * GasConstant * T);
1090}
1091
1092double HMWSoln::ADebye_J(double tempArg, double presArg) const
1093{
1094 double T = temperature();
1095 if (tempArg != -1.0) {
1096 T = tempArg;
1097 }
1098 double A_L = ADebye_L(T, presArg);
1099 double d2 = d2A_DebyedT2_TP(T, presArg);
1100 double d2Aphi = d2 / 3.0;
1101 return 2.0 * A_L / T + 4.0 * GasConstant * T * T *d2Aphi;
1102}
1103
1104double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
1105{
1106 double T = temperature();
1107 if (tempArg != -1.0) {
1108 T = tempArg;
1109 }
1110 double P = pressure();
1111 if (presArg != -1.0) {
1112 P = presArg;
1113 }
1114 double d2AdT2;
1115 switch (m_form_A_Debye) {
1116 case A_DEBYE_CONST:
1117 d2AdT2 = 0.0;
1118 break;
1119 case A_DEBYE_WATER:
1120 d2AdT2 = m_waterProps->ADebye(T, P, 2);
1121 break;
1122 default:
1123 throw CanteraError("HMWSoln::d2A_DebyedT2_TP", "shouldn't be here");
1124 }
1125 return d2AdT2;
1126}
1127
1128// ---------- Other Property Functions
1129
1130// ------------ Private and Restricted Functions ------------------
1131
1133{
1134 m_tmpV.resize(m_kk, 0.0);
1135 m_molalitiesCropped.resize(m_kk, 0.0);
1136
1137 size_t maxCounterIJlen = 1 + (m_kk-1) * (m_kk-2) / 2;
1138
1139 // Figure out the size of the temperature coefficient arrays
1140 int TCoeffLength = 1;
1141 if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1142 TCoeffLength = 2;
1143 } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1144 TCoeffLength = 5;
1145 }
1146
1147 m_Beta0MX_ij.resize(maxCounterIJlen, 0.0);
1148 m_Beta0MX_ij_L.resize(maxCounterIJlen, 0.0);
1149 m_Beta0MX_ij_LL.resize(maxCounterIJlen, 0.0);
1150 m_Beta0MX_ij_P.resize(maxCounterIJlen, 0.0);
1151 m_Beta0MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1152
1153 m_Beta1MX_ij.resize(maxCounterIJlen, 0.0);
1154 m_Beta1MX_ij_L.resize(maxCounterIJlen, 0.0);
1155 m_Beta1MX_ij_LL.resize(maxCounterIJlen, 0.0);
1156 m_Beta1MX_ij_P.resize(maxCounterIJlen, 0.0);
1157 m_Beta1MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1158
1159 m_Beta2MX_ij.resize(maxCounterIJlen, 0.0);
1160 m_Beta2MX_ij_L.resize(maxCounterIJlen, 0.0);
1161 m_Beta2MX_ij_LL.resize(maxCounterIJlen, 0.0);
1162 m_Beta2MX_ij_P.resize(maxCounterIJlen, 0.0);
1163 m_Beta2MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1164
1165 m_CphiMX_ij.resize(maxCounterIJlen, 0.0);
1166 m_CphiMX_ij_L.resize(maxCounterIJlen, 0.0);
1167 m_CphiMX_ij_LL.resize(maxCounterIJlen, 0.0);
1168 m_CphiMX_ij_P.resize(maxCounterIJlen, 0.0);
1169 m_CphiMX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1170
1171 m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
1172 m_Alpha2MX_ij.resize(maxCounterIJlen, 12.0);
1173 m_Theta_ij.resize(maxCounterIJlen, 0.0);
1174 m_Theta_ij_L.resize(maxCounterIJlen, 0.0);
1175 m_Theta_ij_LL.resize(maxCounterIJlen, 0.0);
1176 m_Theta_ij_P.resize(maxCounterIJlen, 0.0);
1177 m_Theta_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1178
1179 size_t n = m_kk*m_kk*m_kk;
1180 m_Psi_ijk.resize(n, 0.0);
1181 m_Psi_ijk_L.resize(n, 0.0);
1182 m_Psi_ijk_LL.resize(n, 0.0);
1183 m_Psi_ijk_P.resize(n, 0.0);
1184 m_Psi_ijk_coeff.resize(TCoeffLength, n, 0.0);
1185
1186 m_Lambda_nj.resize(m_kk, m_kk, 0.0);
1190 m_Lambda_nj_coeff.resize(TCoeffLength, m_kk * m_kk, 0.0);
1191
1192 m_Mu_nnn.resize(m_kk, 0.0);
1193 m_Mu_nnn_L.resize(m_kk, 0.0);
1194 m_Mu_nnn_LL.resize(m_kk, 0.0);
1195 m_Mu_nnn_P.resize(m_kk, 0.0);
1196 m_Mu_nnn_coeff.resize(TCoeffLength, m_kk, 0.0);
1197
1198 m_lnActCoeffMolal_Scaled.resize(m_kk, 0.0);
1199 m_dlnActCoeffMolaldT_Scaled.resize(m_kk, 0.0);
1201 m_dlnActCoeffMolaldP_Scaled.resize(m_kk, 0.0);
1202 m_lnActCoeffMolal_Unscaled.resize(m_kk, 0.0);
1206
1207 m_CounterIJ.resize(m_kk*m_kk, 0);
1208 m_gfunc_IJ.resize(maxCounterIJlen, 0.0);
1209 m_g2func_IJ.resize(maxCounterIJlen, 0.0);
1210 m_hfunc_IJ.resize(maxCounterIJlen, 0.0);
1211 m_h2func_IJ.resize(maxCounterIJlen, 0.0);
1212 m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1213 m_BMX_IJ_L.resize(maxCounterIJlen, 0.0);
1214 m_BMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1215 m_BMX_IJ_P.resize(maxCounterIJlen, 0.0);
1216 m_BprimeMX_IJ.resize(maxCounterIJlen, 0.0);
1217 m_BprimeMX_IJ_L.resize(maxCounterIJlen, 0.0);
1218 m_BprimeMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1219 m_BprimeMX_IJ_P.resize(maxCounterIJlen, 0.0);
1220 m_BphiMX_IJ.resize(maxCounterIJlen, 0.0);
1221 m_BphiMX_IJ_L.resize(maxCounterIJlen, 0.0);
1222 m_BphiMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1223 m_BphiMX_IJ_P.resize(maxCounterIJlen, 0.0);
1224 m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1225 m_Phi_IJ_L.resize(maxCounterIJlen, 0.0);
1226 m_Phi_IJ_LL.resize(maxCounterIJlen, 0.0);
1227 m_Phi_IJ_P.resize(maxCounterIJlen, 0.0);
1228 m_Phiprime_IJ.resize(maxCounterIJlen, 0.0);
1229 m_PhiPhi_IJ.resize(maxCounterIJlen, 0.0);
1230 m_PhiPhi_IJ_L.resize(maxCounterIJlen, 0.0);
1231 m_PhiPhi_IJ_LL.resize(maxCounterIJlen, 0.0);
1232 m_PhiPhi_IJ_P.resize(maxCounterIJlen, 0.0);
1233 m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1234 m_CMX_IJ_L.resize(maxCounterIJlen, 0.0);
1235 m_CMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1236 m_CMX_IJ_P.resize(maxCounterIJlen, 0.0);
1237
1238 m_gamma_tmp.resize(m_kk, 0.0);
1239 IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
1240 CROP_speciesCropped_.resize(m_kk, 0);
1241
1243}
1244
1246{
1247 static const int cacheId = m_cache.getId();
1248 CachedScalar cached = m_cache.getScalar(cacheId);
1249 if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
1250 return;
1251 }
1252
1253 // Calculate the molalities. Currently, the molalities may not be current
1254 // with respect to the contents of the State objects' data.
1256
1257 // Calculate a cropped set of molalities that will be used in all activity
1258 // coefficient calculations.
1260
1261 // Update the temperature dependence of the Pitzer coefficients and their
1262 // derivatives
1264
1265 // Calculate the IMS cutoff factors
1267
1268 // Now do the main calculation.
1270
1271 double xmolSolvent = moleFraction(0);
1272 double xx = std::max(m_xmolSolventMIN, xmolSolvent);
1273 double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1274 double lnxs = log(xx);
1275
1276 for (size_t k = 1; k < m_kk; k++) {
1277 CROP_speciesCropped_[k] = 0;
1279 if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max- 2.5 *lnxs)) {
1280 CROP_speciesCropped_[k] = 2;
1281 m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max - 2.5 * lnxs;
1282 }
1283 if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
1284 // -1.0 and -1.5 caused multiple solutions
1285 CROP_speciesCropped_[k] = 2;
1286 m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
1287 }
1288 }
1289 CROP_speciesCropped_[0] = 0;
1290 m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
1291 if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
1292 CROP_speciesCropped_[0] = 2;
1293 m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
1294 }
1295 if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
1296 CROP_speciesCropped_[0] = 2;
1297 // -0.5 caused multiple solutions
1298 m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
1299 }
1300 if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max - 0.5 * lnxs) {
1301 CROP_speciesCropped_[0] = 2;
1302 m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max - 0.5 * lnxs;
1303 }
1304
1305 // Now do the pH Scaling
1307}
1308
1310{
1311 double Imax = 0.0;
1312 m_molalitiesAreCropped = false;
1313
1314 for (size_t k = 0; k < m_kk; k++) {
1316 Imax = std::max(m_molalities[k] * charge(k) * charge(k), Imax);
1317 }
1318
1319 int cropMethod = 1;
1320 if (cropMethod == 0) {
1321 // Quick return
1322 if (Imax < m_maxIionicStrength) {
1323 return;
1324 }
1325
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) {
1331 continue;
1332 }
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);
1336
1337 // Only loop over oppositely charge species
1338 if (charge_i * charge_j < 0) {
1339 double Iac_max = m_maxIionicStrength;
1340
1342 Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
1343 if (Imax > Iac_max) {
1344 m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
1345 }
1346 Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
1347 if (Imax > Iac_max) {
1348 m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
1349 }
1350 } else {
1351 Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
1352 if (Imax > Iac_max) {
1353 m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
1354 }
1355 Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
1356 if (Imax > Iac_max) {
1357 m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
1358 }
1359 }
1360 }
1361 }
1362 }
1363
1364 // Do this loop 10 times until we have achieved charge neutrality in the
1365 // cropped molalities
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) {
1376 double anion_contrib = - m_molalitiesCropped[i] * charge_i;
1377 anion_charge += anion_contrib;
1378 if (anion_contrib > anion_contrib_max) {
1379 anion_contrib_max = anion_contrib;
1380 anion_contrib_max_i = i;
1381 }
1382 } else if (charge_i > 0.0) {
1383 double cation_contrib = m_molalitiesCropped[i] * charge_i;
1384 cation_charge += cation_contrib;
1385 if (cation_contrib > cation_contrib_max) {
1386 cation_contrib_max = cation_contrib;
1387 cation_contrib_max_i = i;
1388 }
1389 }
1390 }
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);
1394 double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
1395 if (desiredCrop < maxCrop) {
1396 m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
1397 break;
1398 } else {
1399 m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
1400 }
1401 } else if (total_charge < -1.0E-8) {
1402 double desiredCrop = total_charge/charge(anion_contrib_max_i);
1403 double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
1404 if (desiredCrop < maxCrop) {
1405 m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
1406 break;
1407 } else {
1408 m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
1409 }
1410 } else {
1411 break;
1412 }
1413 }
1414 }
1415
1416 if (cropMethod == 1) {
1417 double* molF = m_gamma_tmp.data();
1418 getMoleFractions(molF);
1419 double xmolSolvent = molF[0];
1420 if (xmolSolvent >= MC_X_o_cutoff_) {
1421 return;
1422 }
1423
1425 double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1426 double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1427 double denomInv = 1.0/ (m_Mnaught * p);
1428 for (size_t k = 0; k < m_kk; k++) {
1429 m_molalitiesCropped[k] = molF[k] * denomInv;
1430 }
1431
1432 // Do a further check to see if the Ionic strength is below a max value
1433 // Reduce the molalities to enforce this. Note, this algorithm preserves
1434 // the charge neutrality of the solution after cropping.
1435 double Itmp = 0.0;
1436 for (size_t k = 0; k < m_kk; k++) {
1437 Itmp += m_molalitiesCropped[k] * charge(k) * charge(k);
1438 }
1439 if (Itmp > m_maxIionicStrength) {
1440 double ratio = Itmp / m_maxIionicStrength;
1441 for (size_t k = 0; k < m_kk; k++) {
1442 if (charge(k) != 0.0) {
1443 m_molalitiesCropped[k] *= ratio;
1444 }
1445 }
1446 }
1447 }
1448}
1449
1451{
1452 m_CounterIJ.resize(m_kk * m_kk);
1453 int counter = 0;
1454 for (size_t i = 0; i < m_kk; i++) {
1455 size_t n = i;
1456 size_t nc = m_kk * i;
1457 m_CounterIJ[n] = 0;
1458 m_CounterIJ[nc] = 0;
1459 }
1460 for (size_t i = 1; i < (m_kk - 1); i++) {
1461 size_t n = m_kk * i + i;
1462 m_CounterIJ[n] = 0;
1463 for (size_t j = (i+1); j < m_kk; j++) {
1464 n = m_kk * j + i;
1465 size_t nc = m_kk * i + j;
1466 counter++;
1467 m_CounterIJ[n] = counter;
1468 m_CounterIJ[nc] = counter;
1469 }
1470 }
1471}
1472
1474{
1475 double IMS_gamma_o_min_ = 1.0E-5; // value at the zero solvent point
1476 double IMS_gamma_k_min_ = 10.0; // minimum at the zero solvent point
1477 double IMS_slopefCut_ = 0.6; // slope of the f function at the zero solvent point
1478
1479 IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
1480 IMS_efCut_ = 0.0;
1481 bool converged = false;
1482 double oldV = 0.0;
1483 for (int its = 0; its < 100 && !converged; its++) {
1484 oldV = IMS_efCut_;
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;
1487 IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
1488 /
1490 double tmp = IMS_afCut_ + IMS_X_o_cutoff_*(IMS_bfCut_ + IMS_dfCut_ *IMS_X_o_cutoff_);
1491 double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1492 IMS_efCut_ = - eterm * tmp;
1493 if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
1494 converged = true;
1495 }
1496 }
1497 if (!converged) {
1498 throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1499 " failed to converge on the f polynomial");
1500 }
1501 converged = false;
1502 double f_0 = IMS_afCut_ + IMS_efCut_;
1503 double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
1504 IMS_egCut_ = 0.0;
1505 for (int its = 0; its < 100 && !converged; its++) {
1506 oldV = IMS_egCut_;
1507 double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
1508 IMS_agCut_ = exp(lng_0) - IMS_egCut_;
1509 IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
1510 IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
1511 /
1513 double tmp = IMS_agCut_ + IMS_X_o_cutoff_*(IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
1514 double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1515 IMS_egCut_ = - eterm * tmp;
1516 if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
1517 converged = true;
1518 }
1519 }
1520 if (!converged) {
1521 throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1522 " failed to converge on the g polynomial");
1523 }
1524}
1525
1527{
1528 double MC_X_o_min_ = 0.35; // value at the zero solvent point
1529 MC_X_o_cutoff_ = 0.6;
1530 double MC_slopepCut_ = 0.02; // slope of the p function at the zero solvent point
1531 MC_cpCut_ = 0.25;
1532
1533 // Initial starting values
1534 MC_apCut_ = MC_X_o_min_;
1535 MC_epCut_ = 0.0;
1536 bool converged = false;
1537 double oldV = 0.0;
1538 double damp = 0.5;
1539 for (int its = 0; its < 500 && !converged; its++) {
1540 oldV = MC_epCut_;
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_)
1545 /
1546 (MC_X_o_cutoff_ * MC_X_o_cutoff_/MC_cpCut_ - 2.0 * MC_X_o_cutoff_));
1547 MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
1548 double tmp = MC_apCut_ + MC_X_o_cutoff_*(MC_bpCut_ + MC_dpCut_ * MC_X_o_cutoff_);
1549 double eterm = std::exp(- MC_X_o_cutoff_ / MC_cpCut_);
1550 MC_epCut_ = - eterm * tmp;
1551 double diff = MC_epCut_ - oldV;
1552 if (fabs(diff) < 1.0E-14) {
1553 converged = true;
1554 }
1555 }
1556 if (!converged) {
1557 throw CanteraError("HMWSoln::calcMCCutoffParams_()",
1558 " failed to converge on the p polynomial");
1559 }
1560}
1561
1563{
1564 double T = temperature();
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;
1570 if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1571 tlin = T - m_TempPitzerRef;
1572 } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1573 tlin = T - m_TempPitzerRef;
1574 tquad = T * T - m_TempPitzerRef * m_TempPitzerRef;
1575 tln = log(T/ m_TempPitzerRef);
1576 tinv = 1.0/T - 1.0/m_TempPitzerRef;
1577 }
1578
1579 for (size_t i = 1; i < (m_kk - 1); i++) {
1580 for (size_t j = (i+1); j < m_kk; j++) {
1581 // Find the counterIJ for the symmetric binary interaction
1582 size_t n = m_kk*i + j;
1583 size_t counterIJ = m_CounterIJ[n];
1584
1585 const double* beta0MX_coeff = m_Beta0MX_ij_coeff.ptrColumn(counterIJ);
1586 const double* beta1MX_coeff = m_Beta1MX_ij_coeff.ptrColumn(counterIJ);
1587 const double* beta2MX_coeff = m_Beta2MX_ij_coeff.ptrColumn(counterIJ);
1588 const double* CphiMX_coeff = m_CphiMX_ij_coeff.ptrColumn(counterIJ);
1589 const double* Theta_coeff = m_Theta_ij_coeff.ptrColumn(counterIJ);
1590
1591 switch (m_formPitzerTemp) {
1592 case PITZER_TEMP_CONSTANT:
1593 break;
1594 case PITZER_TEMP_LINEAR:
1595
1596 m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1597 + beta0MX_coeff[1]*tlin;
1598 m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1];
1599 m_Beta0MX_ij_LL[counterIJ] = 0.0;
1600 m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1601 + beta1MX_coeff[1]*tlin;
1602 m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1];
1603 m_Beta1MX_ij_LL[counterIJ] = 0.0;
1604 m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1605 + beta2MX_coeff[1]*tlin;
1606 m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1];
1607 m_Beta2MX_ij_LL[counterIJ] = 0.0;
1608 m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1609 + CphiMX_coeff[1]*tlin;
1610 m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1];
1611 m_CphiMX_ij_LL[counterIJ] = 0.0;
1612 m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1613 m_Theta_ij_L[counterIJ] = Theta_coeff[1];
1614 m_Theta_ij_LL[counterIJ] = 0.0;
1615 break;
1616
1617 case PITZER_TEMP_COMPLEX1:
1618 m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1619 + beta0MX_coeff[1]*tlin
1620 + beta0MX_coeff[2]*tquad
1621 + beta0MX_coeff[3]*tinv
1622 + beta0MX_coeff[4]*tln;
1623 m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1624 + beta1MX_coeff[1]*tlin
1625 + beta1MX_coeff[2]*tquad
1626 + beta1MX_coeff[3]*tinv
1627 + beta1MX_coeff[4]*tln;
1628 m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1629 + beta2MX_coeff[1]*tlin
1630 + beta2MX_coeff[2]*tquad
1631 + beta2MX_coeff[3]*tinv
1632 + beta2MX_coeff[4]*tln;
1633 m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1634 + CphiMX_coeff[1]*tlin
1635 + CphiMX_coeff[2]*tquad
1636 + CphiMX_coeff[3]*tinv
1637 + CphiMX_coeff[4]*tln;
1638 m_Theta_ij[counterIJ] = Theta_coeff[0]
1639 + Theta_coeff[1]*tlin
1640 + Theta_coeff[2]*tquad
1641 + Theta_coeff[3]*tinv
1642 + Theta_coeff[4]*tln;
1643 m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1]
1644 + beta0MX_coeff[2]*twoT
1645 - beta0MX_coeff[3]*invT2
1646 + beta0MX_coeff[4]*invT;
1647 m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1]
1648 + beta1MX_coeff[2]*twoT
1649 - beta1MX_coeff[3]*invT2
1650 + beta1MX_coeff[4]*invT;
1651 m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1]
1652 + beta2MX_coeff[2]*twoT
1653 - beta2MX_coeff[3]*invT2
1654 + beta2MX_coeff[4]*invT;
1655 m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1]
1656 + CphiMX_coeff[2]*twoT
1657 - CphiMX_coeff[3]*invT2
1658 + CphiMX_coeff[4]*invT;
1659 m_Theta_ij_L[counterIJ] = Theta_coeff[1]
1660 + Theta_coeff[2]*twoT
1661 - Theta_coeff[3]*invT2
1662 + Theta_coeff[4]*invT;
1663 doDerivs = 2;
1664 if (doDerivs > 1) {
1665 m_Beta0MX_ij_LL[counterIJ] =
1666 + beta0MX_coeff[2]*2.0
1667 + beta0MX_coeff[3]*twoinvT3
1668 - beta0MX_coeff[4]*invT2;
1669 m_Beta1MX_ij_LL[counterIJ] =
1670 + beta1MX_coeff[2]*2.0
1671 + beta1MX_coeff[3]*twoinvT3
1672 - beta1MX_coeff[4]*invT2;
1673 m_Beta2MX_ij_LL[counterIJ] =
1674 + beta2MX_coeff[2]*2.0
1675 + beta2MX_coeff[3]*twoinvT3
1676 - beta2MX_coeff[4]*invT2;
1677 m_CphiMX_ij_LL[counterIJ] =
1678 + CphiMX_coeff[2]*2.0
1679 + CphiMX_coeff[3]*twoinvT3
1680 - CphiMX_coeff[4]*invT2;
1681 m_Theta_ij_LL[counterIJ] =
1682 + Theta_coeff[2]*2.0
1683 + Theta_coeff[3]*twoinvT3
1684 - Theta_coeff[4]*invT2;
1685 }
1686 break;
1687 }
1688 }
1689 }
1690
1691 // Lambda interactions and Mu_nnn
1692 // i must be neutral for this term to be nonzero. We take advantage of this
1693 // here to lower the operation count.
1694 for (size_t i = 1; i < m_kk; i++) {
1695 if (charge(i) == 0.0) {
1696 for (size_t j = 1; j < m_kk; j++) {
1697 size_t n = i * m_kk + j;
1698 const double* Lambda_coeff = m_Lambda_nj_coeff.ptrColumn(n);
1699 switch (m_formPitzerTemp) {
1700 case PITZER_TEMP_CONSTANT:
1701 m_Lambda_nj(i,j) = Lambda_coeff[0];
1702 break;
1703 case PITZER_TEMP_LINEAR:
1704 m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1705 m_Lambda_nj_L(i,j) = Lambda_coeff[1];
1706 m_Lambda_nj_LL(i,j) = 0.0;
1707 break;
1708 case PITZER_TEMP_COMPLEX1:
1709 m_Lambda_nj(i,j) = Lambda_coeff[0]
1710 + Lambda_coeff[1]*tlin
1711 + Lambda_coeff[2]*tquad
1712 + Lambda_coeff[3]*tinv
1713 + Lambda_coeff[4]*tln;
1714
1715 m_Lambda_nj_L(i,j) = Lambda_coeff[1]
1716 + Lambda_coeff[2]*twoT
1717 - Lambda_coeff[3]*invT2
1718 + Lambda_coeff[4]*invT;
1719
1720 m_Lambda_nj_LL(i,j) =
1721 Lambda_coeff[2]*2.0
1722 + Lambda_coeff[3]*twoinvT3
1723 - Lambda_coeff[4]*invT2;
1724 }
1725
1726 if (j == i) {
1727 const double* Mu_coeff = m_Mu_nnn_coeff.ptrColumn(i);
1728 switch (m_formPitzerTemp) {
1729 case PITZER_TEMP_CONSTANT:
1730 m_Mu_nnn[i] = Mu_coeff[0];
1731 break;
1732 case PITZER_TEMP_LINEAR:
1733 m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1734 m_Mu_nnn_L[i] = Mu_coeff[1];
1735 m_Mu_nnn_LL[i] = 0.0;
1736 break;
1737 case PITZER_TEMP_COMPLEX1:
1738 m_Mu_nnn[i] = Mu_coeff[0]
1739 + Mu_coeff[1]*tlin
1740 + Mu_coeff[2]*tquad
1741 + Mu_coeff[3]*tinv
1742 + Mu_coeff[4]*tln;
1743 m_Mu_nnn_L[i] = Mu_coeff[1]
1744 + Mu_coeff[2]*twoT
1745 - Mu_coeff[3]*invT2
1746 + Mu_coeff[4]*invT;
1747 m_Mu_nnn_LL[i] =
1748 Mu_coeff[2]*2.0
1749 + Mu_coeff[3]*twoinvT3
1750 - Mu_coeff[4]*invT2;
1751 }
1752 }
1753 }
1754 }
1755 }
1756
1757 switch(m_formPitzerTemp) {
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++) {
1762 size_t n = i * m_kk *m_kk + j * m_kk + k;
1763 const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1764 m_Psi_ijk[n] = Psi_coeff[0];
1765 }
1766 }
1767 }
1768 break;
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++) {
1773 size_t n = i * m_kk *m_kk + j * m_kk + k;
1774 const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1775 m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1776 m_Psi_ijk_L[n] = Psi_coeff[1];
1777 m_Psi_ijk_LL[n] = 0.0;
1778 }
1779 }
1780 }
1781 break;
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++) {
1786 size_t n = i * m_kk *m_kk + j * m_kk + k;
1787 const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1788 m_Psi_ijk[n] = Psi_coeff[0]
1789 + Psi_coeff[1]*tlin
1790 + Psi_coeff[2]*tquad
1791 + Psi_coeff[3]*tinv
1792 + Psi_coeff[4]*tln;
1793 m_Psi_ijk_L[n] = Psi_coeff[1]
1794 + Psi_coeff[2]*twoT
1795 - Psi_coeff[3]*invT2
1796 + Psi_coeff[4]*invT;
1797 m_Psi_ijk_LL[n] =
1798 Psi_coeff[2]*2.0
1799 + Psi_coeff[3]*twoinvT3
1800 - Psi_coeff[4]*invT2;
1801 }
1802 }
1803 }
1804 break;
1805 }
1806}
1807
1809{
1810 // Use the CROPPED molality of the species in solution.
1811 const vector<double>& molality = m_molalitiesCropped;
1812
1813 // These are data inputs about the Pitzer correlation. They come from the
1814 // input file for the Pitzer model.
1815 vector<double>& gamma_Unscaled = m_gamma_tmp;
1816
1817 // Local variables defined by Coltrin
1818 double etheta[5][5], etheta_prime[5][5], sqrtIs;
1819
1820 // Molality based ionic strength of the solution
1821 double Is = 0.0;
1822
1823 // Molar charge of the solution: In Pitzer's notation, this is his variable
1824 // called "Z".
1825 double molarcharge = 0.0;
1826
1827 // molalitysum is the sum of the molalities over all solutes, even those
1828 // with zero charge.
1829 double molalitysumUncropped = 0.0;
1830
1831 // Make sure the counter variables are setup
1833
1834 // ---------- Calculate common sums over solutes ---------------------
1835 for (size_t n = 1; n < m_kk; n++) {
1836 // ionic strength
1837 Is += charge(n) * charge(n) * molality[n];
1838 // total molar charge
1839 molarcharge += fabs(charge(n)) * molality[n];
1840 molalitysumUncropped += m_molalities[n];
1841 }
1842 Is *= 0.5;
1843
1844 // Store the ionic molality in the object for reference.
1845 m_IionicMolality = Is;
1846 sqrtIs = sqrt(Is);
1847
1848 // The following call to calc_lambdas() calculates all 16 elements of the
1849 // elambda and elambda1 arrays, given the value of the ionic strength (Is)
1850 calc_lambdas(Is);
1851
1852 // Step 2: Find the coefficients E-theta and E-thetaprime for all
1853 // combinations of positive unlike charges up to 4
1854 for (int z1 = 1; z1 <=4; z1++) {
1855 for (int z2 =1; z2 <=4; z2++) {
1856 calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
1857 }
1858 }
1859
1860 // calculate g(x) and hfunc(x) for each cation-anion pair MX. In the
1861 // original literature, hfunc, was called gprime. However, it's not the
1862 // derivative of g(x), so I renamed it.
1863 for (size_t i = 1; i < (m_kk - 1); i++) {
1864 for (size_t j = (i+1); j < m_kk; j++) {
1865 // Find the counterIJ for the symmetric binary interaction
1866 size_t n = m_kk*i + j;
1867 size_t counterIJ = m_CounterIJ[n];
1868
1869 // Only loop over oppositely charge species
1870 if (charge(i)*charge(j) < 0) {
1871 // x is a reduced function variable
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);
1875 m_hfunc_IJ[counterIJ] = -2.0 *
1876 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
1877 } else {
1878 m_gfunc_IJ[counterIJ] = 0.0;
1879 m_hfunc_IJ[counterIJ] = 0.0;
1880 }
1881
1882 if (m_Beta2MX_ij[counterIJ] != 0.0) {
1883 double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
1884 if (x2 > 1.0E-100) {
1885 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
1886 m_h2func_IJ[counterIJ] = -2.0 *
1887 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
1888 } else {
1889 m_g2func_IJ[counterIJ] = 0.0;
1890 m_h2func_IJ[counterIJ] = 0.0;
1891 }
1892 }
1893 } else {
1894 m_gfunc_IJ[counterIJ] = 0.0;
1895 m_hfunc_IJ[counterIJ] = 0.0;
1896 }
1897 }
1898 }
1899
1900 // SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX
1901 // Agrees with Pitzer, Eq. (49), (51), (55)
1902 for (size_t i = 1; i < m_kk - 1; i++) {
1903 for (size_t j = i+1; j < m_kk; j++) {
1904 // Find the counterIJ for the symmetric binary interaction
1905 size_t n = m_kk*i + j;
1906 size_t counterIJ = m_CounterIJ[n];
1907
1908 // both species have a non-zero charge, and one is positive and the
1909 // other is negative
1910 if (charge(i)*charge(j) < 0.0) {
1911 m_BMX_IJ[counterIJ] = m_Beta0MX_ij[counterIJ]
1912 + m_Beta1MX_ij[counterIJ] * m_gfunc_IJ[counterIJ]
1913 + m_Beta2MX_ij[counterIJ] * m_g2func_IJ[counterIJ];
1914
1915 if (Is > 1.0E-150) {
1916 m_BprimeMX_IJ[counterIJ] = (m_Beta1MX_ij[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
1917 m_Beta2MX_ij[counterIJ] * m_h2func_IJ[counterIJ]/Is);
1918 } else {
1919 m_BprimeMX_IJ[counterIJ] = 0.0;
1920 }
1921 m_BphiMX_IJ[counterIJ] = m_BMX_IJ[counterIJ] + Is*m_BprimeMX_IJ[counterIJ];
1922 } else {
1923 m_BMX_IJ[counterIJ] = 0.0;
1924 m_BprimeMX_IJ[counterIJ] = 0.0;
1925 m_BphiMX_IJ[counterIJ] = 0.0;
1926 }
1927 }
1928 }
1929
1930 // SUBSECTION TO CALCULATE CMX
1931 // Agrees with Pitzer, Eq. (53).
1932 for (size_t i = 1; i < m_kk-1; i++) {
1933 for (size_t j = i+1; j < m_kk; j++) {
1934 // Find the counterIJ for the symmetric binary interaction
1935 size_t n = m_kk*i + j;
1936 size_t counterIJ = m_CounterIJ[n];
1937
1938 // both species have a non-zero charge, and one is positive
1939 // and the other is negative
1940 if (charge(i)*charge(j) < 0.0) {
1941 m_CMX_IJ[counterIJ] = m_CphiMX_ij[counterIJ]/
1942 (2.0* sqrt(fabs(charge(i)*charge(j))));
1943 } else {
1944 m_CMX_IJ[counterIJ] = 0.0;
1945 }
1946 }
1947 }
1948
1949 // SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi
1950 // Agrees with Pitzer, Eq. 72, 73, 74
1951 for (size_t i = 1; i < m_kk-1; i++) {
1952 for (size_t j = i+1; j < m_kk; j++) {
1953 // Find the counterIJ for the symmetric binary interaction
1954 size_t n = m_kk*i + j;
1955 size_t counterIJ = m_CounterIJ[n];
1956
1957 // both species have a non-zero charge, and one is positive and the
1958 // other is negative
1959 if (charge(i)*charge(j) > 0) {
1960 int z1 = (int) fabs(charge(i));
1961 int z2 = (int) fabs(charge(j));
1962 m_Phi_IJ[counterIJ] = m_Theta_ij[counterIJ] + etheta[z1][z2];
1963 m_Phiprime_IJ[counterIJ] = etheta_prime[z1][z2];
1964 m_PhiPhi_IJ[counterIJ] = m_Phi_IJ[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
1965 } else {
1966 m_Phi_IJ[counterIJ] = 0.0;
1967 m_Phiprime_IJ[counterIJ] = 0.0;
1968 m_PhiPhi_IJ[counterIJ] = 0.0;
1969 }
1970 }
1971 }
1972
1973 // SUBSECTION FOR CALCULATION OF F
1974 // Agrees with Pitzer Eqn. (65)
1975 double Aphi = A_Debye_TP() / 3.0;
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++) {
1980 // Find the counterIJ for the symmetric binary interaction
1981 size_t n = m_kk*i + j;
1982 size_t counterIJ = m_CounterIJ[n];
1983
1984 // both species have a non-zero charge, and one is positive and the
1985 // other is negative
1986 if (charge(i)*charge(j) < 0) {
1987 F += molality[i]*molality[j] * m_BprimeMX_IJ[counterIJ];
1988 }
1989
1990 // Both species have a non-zero charge, and they
1991 // have the same sign
1992 if (charge(i)*charge(j) > 0) {
1993 F += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
1994 }
1995 }
1996 }
1997
1998 for (size_t i = 1; i < m_kk; i++) {
1999
2000 // SUBSECTION FOR CALCULATING THE ACTCOEFF FOR CATIONS
2001 // equations agree with my notes, Eqn. (118).
2002 // Equations agree with Pitzer, eqn.(63)
2003 if (charge(i) > 0.0) {
2004 // species i is the cation (positive) to calc the actcoeff
2005 double zsqF = charge(i)*charge(i)*F;
2006 double sum1 = 0.0;
2007 double sum2 = 0.0;
2008 double sum3 = 0.0;
2009 double sum4 = 0.0;
2010 double sum5 = 0.0;
2011 for (size_t j = 1; j < m_kk; j++) {
2012 // Find the counterIJ for the symmetric binary interaction
2013 size_t n = m_kk*i + j;
2014 size_t counterIJ = m_CounterIJ[n];
2015
2016 if (charge(j) < 0.0) {
2017 // sum over all anions
2018 sum1 += molality[j] *
2019 (2.0*m_BMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]);
2020 if (j < m_kk-1) {
2021 // This term is the ternary interaction involving the
2022 // non-duplicate sum over double anions, j, k, with
2023 // respect to the cation, i.
2024 for (size_t k = j+1; k < m_kk; k++) {
2025 // an inner sum over all anions
2026 if (charge(k) < 0.0) {
2027 n = k + j * m_kk + i * m_kk * m_kk;
2028 sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
2029 }
2030 }
2031 }
2032 }
2033
2034 if (charge(j) > 0.0) {
2035 // sum over all cations
2036 if (j != i) {
2037 sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
2038 }
2039 for (size_t k = 1; k < m_kk; k++) {
2040 if (charge(k) < 0.0) {
2041 // two inner sums over anions
2042 n = k + j * m_kk + i * m_kk * m_kk;
2043 sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
2044
2045 // Find the counterIJ for the j,k interaction
2046 n = m_kk*j + k;
2047 size_t counterIJ2 = m_CounterIJ[n];
2048 sum4 += (fabs(charge(i))*
2049 molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
2050 }
2051 }
2052 }
2053
2054 // Handle neutral j species
2055 if (charge(j) == 0) {
2056 sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
2057
2058 // Zeta interaction term
2059 for (size_t k = 1; k < m_kk; k++) {
2060 if (charge(k) < 0.0) {
2061 size_t izeta = j;
2062 size_t jzeta = i;
2063 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2064 double zeta = m_Psi_ijk[n];
2065 if (zeta != 0.0) {
2066 sum5 += molality[j]*molality[k]*zeta;
2067 }
2068 }
2069 }
2070 }
2071 }
2072
2073 // Add all of the contributions up to yield the log of the solute
2074 // activity coefficients (molality scale)
2075 m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2076 gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2077 }
2078
2079 // SUBSECTION FOR CALCULATING THE ACTCOEFF FOR ANIONS
2080 // equations agree with my notes, Eqn. (119).
2081 // Equations agree with Pitzer, eqn.(64)
2082 if (charge(i) < 0) {
2083 // species i is an anion (negative)
2084 double zsqF = charge(i)*charge(i)*F;
2085 double sum1 = 0.0;
2086 double sum2 = 0.0;
2087 double sum3 = 0.0;
2088 double sum4 = 0.0;
2089 double sum5 = 0.0;
2090 for (size_t j = 1; j < m_kk; j++) {
2091 // Find the counterIJ for the symmetric binary interaction
2092 size_t n = m_kk*i + j;
2093 size_t counterIJ = m_CounterIJ[n];
2094
2095 // For Anions, do the cation interactions.
2096 if (charge(j) > 0) {
2097 sum1 += molality[j]*
2098 (2.0*m_BMX_IJ[counterIJ]+molarcharge*m_CMX_IJ[counterIJ]);
2099 if (j < m_kk-1) {
2100 for (size_t k = j+1; k < m_kk; k++) {
2101 // an inner sum over all cations
2102 if (charge(k) > 0) {
2103 n = k + j * m_kk + i * m_kk * m_kk;
2104 sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
2105 }
2106 }
2107 }
2108 }
2109
2110 // For Anions, do the other anion interactions.
2111 if (charge(j) < 0.0) {
2112 // sum over all anions
2113 if (j != i) {
2114 sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
2115 }
2116 for (size_t k = 1; k < m_kk; k++) {
2117 if (charge(k) > 0.0) {
2118 // two inner sums over cations
2119 n = k + j * m_kk + i * m_kk * m_kk;
2120 sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
2121 // Find the counterIJ for the symmetric binary interaction
2122 n = m_kk*j + k;
2123 size_t counterIJ2 = m_CounterIJ[n];
2124 sum4 += fabs(charge(i))*
2125 molality[j]*molality[k]*m_CMX_IJ[counterIJ2];
2126 }
2127 }
2128 }
2129
2130 // for Anions, do the neutral species interaction
2131 if (charge(j) == 0.0) {
2132 sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
2133 // Zeta interaction term
2134 for (size_t k = 1; k < m_kk; k++) {
2135 if (charge(k) > 0.0) {
2136 size_t izeta = j;
2137 size_t jzeta = k;
2138 size_t kzeta = i;
2139 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2140 double zeta = m_Psi_ijk[n];
2141 if (zeta != 0.0) {
2142 sum5 += molality[j]*molality[k]*zeta;
2143 }
2144 }
2145 }
2146 }
2147 }
2148 m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2149 gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2150 }
2151
2152 // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
2153 // equations agree with my notes,
2154 // Equations agree with Pitzer,
2155 if (charge(i) == 0.0) {
2156 double sum1 = 0.0;
2157 double sum3 = 0.0;
2158 for (size_t j = 1; j < m_kk; j++) {
2159 sum1 += molality[j]*2.0*m_Lambda_nj(i,j);
2160 // Zeta term -> we piggyback on the psi term
2161 if (charge(j) > 0.0) {
2162 for (size_t k = 1; k < m_kk; k++) {
2163 if (charge(k) < 0.0) {
2164 size_t n = k + j * m_kk + i * m_kk * m_kk;
2165 sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
2166 }
2167 }
2168 }
2169 }
2170 double sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i];
2171 m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3;
2172 gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2173 }
2174 }
2175
2176 // SUBSECTION FOR CALCULATING THE OSMOTIC COEFF
2177 // equations agree with my notes, Eqn. (117).
2178 // Equations agree with Pitzer, eqn.(62)
2179 double sum1 = 0.0;
2180 double sum2 = 0.0;
2181 double sum3 = 0.0;
2182 double sum4 = 0.0;
2183 double sum5 = 0.0;
2184 double sum6 = 0.0;
2185 double sum7 = 0.0;
2186
2187 // term1 is the DH term in the osmotic coefficient expression
2188 // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
2189 // implementations.
2190 // Is = Ionic strength on the molality scale (units of (gmol/kg))
2191 // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2192 double term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2193
2194 for (size_t j = 1; j < m_kk; j++) {
2195 // Loop Over Cations
2196 if (charge(j) > 0.0) {
2197 for (size_t k = 1; k < m_kk; k++) {
2198 if (charge(k) < 0.0) {
2199 // Find the counterIJ for the symmetric j,k binary interaction
2200 size_t n = m_kk*j + k;
2201 size_t counterIJ = m_CounterIJ[n];
2202
2203 sum1 += molality[j]*molality[k]*
2204 (m_BphiMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]);
2205 }
2206 }
2207
2208 for (size_t k = j+1; k < m_kk; k++) {
2209 if (j == (m_kk-1)) {
2210 // we should never reach this step
2211 throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2212 "logic error 1 in Step 9 of hmw_act");
2213 }
2214 if (charge(k) > 0.0) {
2215 // Find the counterIJ for the symmetric j,k binary interaction
2216 // between 2 cations.
2217 size_t n = m_kk*j + k;
2218 size_t counterIJ = m_CounterIJ[n];
2219 sum2 += molality[j]*molality[k]*m_PhiPhi_IJ[counterIJ];
2220 for (size_t m = 1; m < m_kk; m++) {
2221 if (charge(m) < 0.0) {
2222 // species m is an anion
2223 n = m + k * m_kk + j * m_kk * m_kk;
2224 sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk[n];
2225 }
2226 }
2227 }
2228 }
2229 }
2230
2231 // Loop Over Anions
2232 if (charge(j) < 0) {
2233 for (size_t k = j+1; k < m_kk; k++) {
2234 if (j == m_kk-1) {
2235 // we should never reach this step
2236 throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2237 "logic error 2 in Step 9 of hmw_act");
2238 }
2239 if (charge(k) < 0) {
2240 // Find the counterIJ for the symmetric j,k binary interaction
2241 // between two anions
2242 size_t n = m_kk*j + k;
2243 size_t counterIJ = m_CounterIJ[n];
2244 sum3 += molality[j]*molality[k]*m_PhiPhi_IJ[counterIJ];
2245 for (size_t m = 1; m < m_kk; m++) {
2246 if (charge(m) > 0.0) {
2247 n = m + k * m_kk + j * m_kk * m_kk;
2248 sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk[n];
2249 }
2250 }
2251 }
2252 }
2253 }
2254
2255 // Loop Over Neutral Species
2256 if (charge(j) == 0) {
2257 for (size_t k = 1; k < m_kk; k++) {
2258 if (charge(k) < 0.0) {
2259 sum4 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2260 }
2261 if (charge(k) > 0.0) {
2262 sum5 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2263 }
2264 if (charge(k) == 0.0) {
2265 if (k > j) {
2266 sum6 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2267 } else if (k == j) {
2268 sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj(j,k);
2269 }
2270 }
2271 if (charge(k) < 0.0) {
2272 size_t izeta = j;
2273 for (size_t m = 1; m < m_kk; m++) {
2274 if (charge(m) > 0.0) {
2275 size_t jzeta = m;
2276 size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2277 double zeta = m_Psi_ijk[n];
2278 if (zeta != 0.0) {
2279 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2280 }
2281 }
2282 }
2283 }
2284 }
2285 sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn[j];
2286 }
2287 }
2288 double sum_m_phi_minus_1 = 2.0 *
2289 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2290 // Calculate the osmotic coefficient from
2291 // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
2292 double osmotic_coef;
2293 if (molalitysumUncropped > 1.0E-150) {
2294 osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2295 } else {
2296 osmotic_coef = 1.0;
2297 }
2298 double lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2299
2300 // In Cantera, we define the activity coefficient of the solvent as
2301 //
2302 // act_0 = actcoeff_0 * Xmol_0
2303 //
2304 // We have just computed act_0. However, this routine returns
2305 // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2306 double xmolSolvent = moleFraction(0);
2307 double xx = std::max(m_xmolSolventMIN, xmolSolvent);
2308 m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
2309}
2310
2312{
2313 static const int cacheId = m_cache.getId();
2314 CachedScalar cached = m_cache.getScalar(cacheId);
2315 if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
2316 return;
2317 }
2318
2319 // Zero the unscaled 2nd derivatives
2321
2322 // Do the actual calculation of the unscaled temperature derivatives
2324
2325 for (size_t k = 1; k < m_kk; k++) {
2326 if (CROP_speciesCropped_[k] == 2) {
2328 }
2329 }
2330
2331 if (CROP_speciesCropped_[0]) {
2333 }
2334
2335 // Do the pH scaling to the derivatives
2337}
2338
2340{
2341 // It may be assumed that the Pitzer activity coefficient routine is called
2342 // immediately preceding the calling of this routine. Therefore, some
2343 // quantities do not need to be recalculated in this routine.
2344
2345 const vector<double>& molality = m_molalitiesCropped;
2346 double* d_gamma_dT_Unscaled = m_gamma_tmp.data();
2347
2348 // Local variables defined by Coltrin
2349 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2350
2351 // Molality based ionic strength of the solution
2352 double Is = 0.0;
2353
2354 // Molar charge of the solution: In Pitzer's notation, this is his variable
2355 // called "Z".
2356 double molarcharge = 0.0;
2357
2358 // molalitysum is the sum of the molalities over all solutes, even those
2359 // with zero charge.
2360 double molalitysum = 0.0;
2361
2362 // Make sure the counter variables are setup
2364
2365 // ---------- Calculate common sums over solutes ---------------------
2366 for (size_t n = 1; n < m_kk; n++) {
2367 // ionic strength
2368 Is += charge(n) * charge(n) * molality[n];
2369 // total molar charge
2370 molarcharge += fabs(charge(n)) * molality[n];
2371 molalitysum += molality[n];
2372 }
2373 Is *= 0.5;
2374
2375 // Store the ionic molality in the object for reference.
2376 m_IionicMolality = Is;
2377 sqrtIs = sqrt(Is);
2378
2379 // The following call to calc_lambdas() calculates all 16 elements of the
2380 // elambda and elambda1 arrays, given the value of the ionic strength (Is)
2381 calc_lambdas(Is);
2382
2383 // Step 2: Find the coefficients E-theta and E-thetaprime for all
2384 // combinations of positive unlike charges up to 4
2385 for (int z1 = 1; z1 <=4; z1++) {
2386 for (int z2 =1; z2 <=4; z2++) {
2387 calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2388 }
2389 }
2390
2391 // calculate g(x) and hfunc(x) for each cation-anion pair MX
2392 // In the original literature, hfunc, was called gprime. However,
2393 // it's not the derivative of g(x), so I renamed it.
2394 for (size_t i = 1; i < (m_kk - 1); i++) {
2395 for (size_t j = (i+1); j < m_kk; j++) {
2396 // Find the counterIJ for the symmetric binary interaction
2397 size_t n = m_kk*i + j;
2398 size_t counterIJ = m_CounterIJ[n];
2399
2400 // Only loop over oppositely charge species
2401 if (charge(i)*charge(j) < 0) {
2402 // x is a reduced function variable
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);
2406 m_hfunc_IJ[counterIJ] = -2.0 *
2407 (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
2408 } else {
2409 m_gfunc_IJ[counterIJ] = 0.0;
2410 m_hfunc_IJ[counterIJ] = 0.0;
2411 }
2412
2413 if (m_Beta2MX_ij_L[counterIJ] != 0.0) {
2414 double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
2415 if (x2 > 1.0E-100) {
2416 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2417 m_h2func_IJ[counterIJ] = -2.0 *
2418 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2419 } else {
2420 m_g2func_IJ[counterIJ] = 0.0;
2421 m_h2func_IJ[counterIJ] = 0.0;
2422 }
2423 }
2424 } else {
2425 m_gfunc_IJ[counterIJ] = 0.0;
2426 m_hfunc_IJ[counterIJ] = 0.0;
2427 }
2428 }
2429 }
2430
2431 // SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L
2432 // These are now temperature derivatives of the previously calculated
2433 // quantities.
2434 for (size_t i = 1; i < m_kk - 1; i++) {
2435 for (size_t j = i+1; j < m_kk; j++) {
2436 // Find the counterIJ for the symmetric binary interaction
2437 size_t n = m_kk*i + j;
2438 size_t counterIJ = m_CounterIJ[n];
2439
2440 // both species have a non-zero charge, and one is positive
2441 // and the other is negative
2442 if (charge(i)*charge(j) < 0.0) {
2443 m_BMX_IJ_L[counterIJ] = m_Beta0MX_ij_L[counterIJ]
2444 + m_Beta1MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ]
2445 + m_Beta2MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ];
2446 if (Is > 1.0E-150) {
2447 m_BprimeMX_IJ_L[counterIJ] = (m_Beta1MX_ij_L[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
2448 m_Beta2MX_ij_L[counterIJ] * m_h2func_IJ[counterIJ]/Is);
2449 } else {
2450 m_BprimeMX_IJ_L[counterIJ] = 0.0;
2451 }
2452 m_BphiMX_IJ_L[counterIJ] = m_BMX_IJ_L[counterIJ] + Is*m_BprimeMX_IJ_L[counterIJ];
2453 } else {
2454 m_BMX_IJ_L[counterIJ] = 0.0;
2455 m_BprimeMX_IJ_L[counterIJ] = 0.0;
2456 m_BphiMX_IJ_L[counterIJ] = 0.0;
2457 }
2458 }
2459 }
2460
2461 // --------- SUBSECTION TO CALCULATE CMX_L ----------
2462 for (size_t i = 1; i < m_kk-1; i++) {
2463 for (size_t j = i+1; j < m_kk; j++) {
2464 // Find the counterIJ for the symmetric binary interaction
2465 size_t n = m_kk*i + j;
2466 size_t counterIJ = m_CounterIJ[n];
2467
2468 // both species have a non-zero charge, and one is positive
2469 // and the other is negative
2470 if (charge(i)*charge(j) < 0.0) {
2471 m_CMX_IJ_L[counterIJ] = m_CphiMX_ij_L[counterIJ]/
2472 (2.0* sqrt(fabs(charge(i)*charge(j))));
2473 } else {
2474 m_CMX_IJ_L[counterIJ] = 0.0;
2475 }
2476 }
2477 }
2478
2479 // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2480 for (size_t i = 1; i < m_kk-1; i++) {
2481 for (size_t j = i+1; j < m_kk; j++) {
2482 // Find the counterIJ for the symmetric binary interaction
2483 size_t n = m_kk*i + j;
2484 size_t counterIJ = m_CounterIJ[n];
2485
2486 // both species have a non-zero charge, and one is positive
2487 // and the other is negative
2488 if (charge(i)*charge(j) > 0) {
2489 m_Phi_IJ_L[counterIJ] = m_Theta_ij_L[counterIJ];
2490 m_Phiprime_IJ[counterIJ] = 0.0;
2491 m_PhiPhi_IJ_L[counterIJ] = m_Phi_IJ_L[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
2492 } else {
2493 m_Phi_IJ_L[counterIJ] = 0.0;
2494 m_Phiprime_IJ[counterIJ] = 0.0;
2495 m_PhiPhi_IJ_L[counterIJ] = 0.0;
2496 }
2497 }
2498 }
2499
2500 // ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
2501 double dA_DebyedT = dA_DebyedT_TP();
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++) {
2507 // Find the counterIJ for the symmetric binary interaction
2508 size_t n = m_kk*i + j;
2509 size_t counterIJ = m_CounterIJ[n];
2510
2511 // both species have a non-zero charge, and one is positive
2512 // and the other is negative
2513 if (charge(i)*charge(j) < 0) {
2514 dFdT += molality[i]*molality[j] * m_BprimeMX_IJ_L[counterIJ];
2515 }
2516
2517 // Both species have a non-zero charge, and they
2518 // have the same sign, that is, both positive or both negative.
2519 if (charge(i)*charge(j) > 0) {
2520 dFdT += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
2521 }
2522 }
2523 }
2524
2525 for (size_t i = 1; i < m_kk; i++) {
2526 // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
2527 if (charge(i) > 0) {
2528 // species i is the cation (positive) to calc the actcoeff
2529 double zsqdFdT = charge(i)*charge(i)*dFdT;
2530 double sum1 = 0.0;
2531 double sum2 = 0.0;
2532 double sum3 = 0.0;
2533 double sum4 = 0.0;
2534 double sum5 = 0.0;
2535 for (size_t j = 1; j < m_kk; j++) {
2536 // Find the counterIJ for the symmetric binary interaction
2537 size_t n = m_kk*i + j;
2538 size_t counterIJ = m_CounterIJ[n];
2539
2540 if (charge(j) < 0.0) {
2541 // sum over all anions
2542 sum1 += molality[j]*
2543 (2.0*m_BMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2544 if (j < m_kk-1) {
2545 // This term is the ternary interaction involving the
2546 // non-duplicate sum over double anions, j, k, with
2547 // respect to the cation, i.
2548 for (size_t k = j+1; k < m_kk; k++) {
2549 // an inner sum over all anions
2550 if (charge(k) < 0.0) {
2551 n = k + j * m_kk + i * m_kk * m_kk;
2552 sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2553 }
2554 }
2555 }
2556 }
2557
2558 if (charge(j) > 0.0) {
2559 // sum over all cations
2560 if (j != i) {
2561 sum2 += molality[j]*(2.0*m_Phi_IJ_L[counterIJ]);
2562 }
2563 for (size_t k = 1; k < m_kk; k++) {
2564 if (charge(k) < 0.0) {
2565 // two inner sums over anions
2566 n = k + j * m_kk + i * m_kk * m_kk;
2567 sum2 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2568
2569 // Find the counterIJ for the j,k interaction
2570 n = m_kk*j + k;
2571 size_t counterIJ2 = m_CounterIJ[n];
2572 sum4 += fabs(charge(i))*
2573 molality[j]*molality[k]*m_CMX_IJ_L[counterIJ2];
2574 }
2575 }
2576 }
2577
2578 // Handle neutral j species
2579 if (charge(j) == 0) {
2580 sum5 += molality[j]*2.0*m_Lambda_nj_L(j,i);
2581 }
2582
2583 // Zeta interaction term
2584 for (size_t k = 1; k < m_kk; k++) {
2585 if (charge(k) < 0.0) {
2586 size_t izeta = j;
2587 size_t jzeta = i;
2588 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2589 double zeta_L = m_Psi_ijk_L[n];
2590 if (zeta_L != 0.0) {
2591 sum5 += molality[j]*molality[k]*zeta_L;
2592 }
2593 }
2594 }
2595 }
2596
2597 // Add all of the contributions up to yield the log of the
2598 // solute activity coefficients (molality scale)
2600 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2601 d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2602 }
2603
2604 // ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------
2605 if (charge(i) < 0) {
2606 // species i is an anion (negative)
2607 double zsqdFdT = charge(i)*charge(i)*dFdT;
2608 double sum1 = 0.0;
2609 double sum2 = 0.0;
2610 double sum3 = 0.0;
2611 double sum4 = 0.0;
2612 double sum5 = 0.0;
2613 for (size_t j = 1; j < m_kk; j++) {
2614 // Find the counterIJ for the symmetric binary interaction
2615 size_t n = m_kk*i + j;
2616 size_t counterIJ = m_CounterIJ[n];
2617
2618 // For Anions, do the cation interactions.
2619 if (charge(j) > 0) {
2620 sum1 += molality[j]*
2621 (2.0*m_BMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2622 if (j < m_kk-1) {
2623 for (size_t k = j+1; k < m_kk; k++) {
2624 // an inner sum over all cations
2625 if (charge(k) > 0) {
2626 n = k + j * m_kk + i * m_kk * m_kk;
2627 sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2628 }
2629 }
2630 }
2631 }
2632
2633 // For Anions, do the other anion interactions.
2634 if (charge(j) < 0.0) {
2635 // sum over all anions
2636 if (j != i) {
2637 sum2 += molality[j]*(2.0*m_Phi_IJ_L[counterIJ]);
2638 }
2639 for (size_t k = 1; k < m_kk; k++) {
2640 if (charge(k) > 0.0) {
2641 // two inner sums over cations
2642 n = k + j * m_kk + i * m_kk * m_kk;
2643 sum2 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2644 // Find the counterIJ for the symmetric binary interaction
2645 n = m_kk*j + k;
2646 size_t counterIJ2 = m_CounterIJ[n];
2647 sum4 += fabs(charge(i)) *
2648 molality[j]*molality[k]*m_CMX_IJ_L[counterIJ2];
2649 }
2650 }
2651 }
2652
2653 // for Anions, do the neutral species interaction
2654 if (charge(j) == 0.0) {
2655 sum5 += molality[j]*2.0*m_Lambda_nj_L(j,i);
2656 for (size_t k = 1; k < m_kk; k++) {
2657 if (charge(k) > 0.0) {
2658 size_t izeta = j;
2659 size_t jzeta = k;
2660 size_t kzeta = i;
2661 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2662 double zeta_L = m_Psi_ijk_L[n];
2663 if (zeta_L != 0.0) {
2664 sum5 += molality[j]*molality[k]*zeta_L;
2665 }
2666 }
2667 }
2668 }
2669 }
2671 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2672 d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2673 }
2674
2675 // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
2676 // equations agree with my notes,
2677 // Equations agree with Pitzer,
2678 if (charge(i) == 0.0) {
2679 double sum1 = 0.0;
2680 double sum3 = 0.0;
2681 for (size_t j = 1; j < m_kk; j++) {
2682 sum1 += molality[j]*2.0*m_Lambda_nj_L(i,j);
2683 // Zeta term -> we piggyback on the psi term
2684 if (charge(j) > 0.0) {
2685 for (size_t k = 1; k < m_kk; k++) {
2686 if (charge(k) < 0.0) {
2687 size_t n = k + j * m_kk + i * m_kk * m_kk;
2688 sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2689 }
2690 }
2691 }
2692 }
2693 double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i];
2694 m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3;
2695 d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2696 }
2697 }
2698
2699 // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT ---------
2700 double sum1 = 0.0;
2701 double sum2 = 0.0;
2702 double sum3 = 0.0;
2703 double sum4 = 0.0;
2704 double sum5 = 0.0;
2705 double sum6 = 0.0;
2706 double sum7 = 0.0;
2707
2708 // term1 is the temperature derivative of the DH term in the osmotic
2709 // coefficient expression
2710 // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
2711 // Is = Ionic strength on the molality scale (units of (gmol/kg))
2712 // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2713 double term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
2714
2715 for (size_t j = 1; j < m_kk; j++) {
2716 // Loop Over Cations
2717 if (charge(j) > 0.0) {
2718 for (size_t k = 1; k < m_kk; k++) {
2719 if (charge(k) < 0.0) {
2720 // Find the counterIJ for the symmetric j,k binary interaction
2721 size_t n = m_kk*j + k;
2722 size_t counterIJ = m_CounterIJ[n];
2723 sum1 += molality[j]*molality[k]*
2724 (m_BphiMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2725 }
2726 }
2727
2728 for (size_t k = j+1; k < m_kk; k++) {
2729 if (j == (m_kk-1)) {
2730 // we should never reach this step
2731 throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2732 "logic error 1 in Step 9 of hmw_act");
2733 }
2734 if (charge(k) > 0.0) {
2735 // Find the counterIJ for the symmetric j,k binary interaction
2736 // between 2 cations.
2737 size_t n = m_kk*j + k;
2738 size_t counterIJ = m_CounterIJ[n];
2739 sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_L[counterIJ];
2740 for (size_t m = 1; m < m_kk; m++) {
2741 if (charge(m) < 0.0) {
2742 // species m is an anion
2743 n = m + k * m_kk + j * m_kk * m_kk;
2744 sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_L[n];
2745 }
2746 }
2747 }
2748 }
2749 }
2750
2751 // Loop Over Anions
2752 if (charge(j) < 0) {
2753 for (size_t k = j+1; k < m_kk; k++) {
2754 if (j == m_kk-1) {
2755 // we should never reach this step
2756 throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2757 "logic error 2 in Step 9 of hmw_act");
2758 }
2759 if (charge(k) < 0) {
2760 // Find the counterIJ for the symmetric j,k binary interaction
2761 // between two anions
2762 size_t n = m_kk*j + k;
2763 size_t counterIJ = m_CounterIJ[n];
2764 sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_L[counterIJ];
2765 for (size_t m = 1; m < m_kk; m++) {
2766 if (charge(m) > 0.0) {
2767 n = m + k * m_kk + j * m_kk * m_kk;
2768 sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_L[n];
2769 }
2770 }
2771 }
2772 }
2773 }
2774
2775 // Loop Over Neutral Species
2776 if (charge(j) == 0) {
2777 for (size_t k = 1; k < m_kk; k++) {
2778 if (charge(k) < 0.0) {
2779 sum4 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2780 }
2781 if (charge(k) > 0.0) {
2782 sum5 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2783 }
2784 if (charge(k) == 0.0) {
2785 if (k > j) {
2786 sum6 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2787 } else if (k == j) {
2788 sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2789 }
2790 }
2791 if (charge(k) < 0.0) {
2792 size_t izeta = j;
2793 for (size_t m = 1; m < m_kk; m++) {
2794 if (charge(m) > 0.0) {
2795 size_t jzeta = m;
2796 size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2797 double zeta_L = m_Psi_ijk_L[n];
2798 if (zeta_L != 0.0) {
2799 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
2800 }
2801 }
2802 }
2803 }
2804 }
2805 sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn_L[j];
2806 }
2807 }
2808 double sum_m_phi_minus_1 = 2.0 *
2809 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2810 // Calculate the osmotic coefficient from
2811 // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
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);
2815 } else {
2816 d_osmotic_coef_dT = 0.0;
2817 }
2818
2819 double d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
2820
2821 // In Cantera, we define the activity coefficient of the solvent as
2822 //
2823 // act_0 = actcoeff_0 * Xmol_0
2824 //
2825 // We have just computed act_0. However, this routine returns
2826 // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2827 m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT;
2828}
2829
2831{
2832 static const int cacheId = m_cache.getId();
2833 CachedScalar cached = m_cache.getScalar(cacheId);
2834 if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
2835 return;
2836 }
2837
2838 // Zero the unscaled 2nd derivatives
2840
2841 //! Calculate the unscaled 2nd derivatives
2843
2844 for (size_t k = 1; k < m_kk; k++) {
2845 if (CROP_speciesCropped_[k] == 2) {
2847 }
2848 }
2849
2850 if (CROP_speciesCropped_[0]) {
2852 }
2853
2854 // Scale the 2nd derivatives
2856}
2857
2859{
2860 const double* molality = m_molalitiesCropped.data();
2861
2862 // Local variables defined by Coltrin
2863 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2864
2865 // Molality based ionic strength of the solution
2866 double Is = 0.0;
2867
2868 // Molar charge of the solution: In Pitzer's notation, this is his variable
2869 // called "Z".
2870 double molarcharge = 0.0;
2871
2872 // molalitysum is the sum of the molalities over all solutes, even those
2873 // with zero charge.
2874 double molalitysum = 0.0;
2875
2876 // Make sure the counter variables are setup
2878
2879 // ---------- Calculate common sums over solutes ---------------------
2880 for (size_t n = 1; n < m_kk; n++) {
2881 // ionic strength
2882 Is += charge(n) * charge(n) * molality[n];
2883 // total molar charge
2884 molarcharge += fabs(charge(n)) * molality[n];
2885 molalitysum += molality[n];
2886 }
2887 Is *= 0.5;
2888
2889 // Store the ionic molality in the object for reference.
2890 m_IionicMolality = Is;
2891 sqrtIs = sqrt(Is);
2892
2893 // The following call to calc_lambdas() calculates all 16 elements of the
2894 // elambda and elambda1 arrays, given the value of the ionic strength (Is)
2895 calc_lambdas(Is);
2896
2897 // Step 2: Find the coefficients E-theta and E-thetaprime for all
2898 // combinations of positive unlike charges up to 4
2899 for (int z1 = 1; z1 <=4; z1++) {
2900 for (int z2 =1; z2 <=4; z2++) {
2901 calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2902 }
2903 }
2904
2905 // calculate gfunc(x) and hfunc(x) for each cation-anion pair MX. In the
2906 // original literature, hfunc, was called gprime. However, it's not the
2907 // derivative of gfunc(x), so I renamed it.
2908 for (size_t i = 1; i < (m_kk - 1); i++) {
2909 for (size_t j = (i+1); j < m_kk; j++) {
2910 // Find the counterIJ for the symmetric binary interaction
2911 size_t n = m_kk*i + j;
2912 size_t counterIJ = m_CounterIJ[n];
2913
2914 // Only loop over oppositely charge species
2915 if (charge(i)*charge(j) < 0) {
2916 // x is a reduced function variable
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);
2920 m_hfunc_IJ[counterIJ] = -2.0*
2921 (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
2922 } else {
2923 m_gfunc_IJ[counterIJ] = 0.0;
2924 m_hfunc_IJ[counterIJ] = 0.0;
2925 }
2926
2927 if (m_Beta2MX_ij_LL[counterIJ] != 0.0) {
2928 double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
2929 if (x2 > 1.0E-100) {
2930 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2931 m_h2func_IJ[counterIJ] = -2.0 *
2932 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2933 } else {
2934 m_g2func_IJ[counterIJ] = 0.0;
2935 m_h2func_IJ[counterIJ] = 0.0;
2936 }
2937 }
2938 } else {
2939 m_gfunc_IJ[counterIJ] = 0.0;
2940 m_hfunc_IJ[counterIJ] = 0.0;
2941 }
2942 }
2943 }
2944
2945 // SUBSECTION TO CALCULATE BMX_L, BprimeMX_LL, BphiMX_L
2946 // These are now temperature derivatives of the previously calculated
2947 // quantities.
2948 for (size_t i = 1; i < m_kk - 1; i++) {
2949 for (size_t j = i+1; j < m_kk; j++) {
2950 // Find the counterIJ for the symmetric binary interaction
2951 size_t n = m_kk*i + j;
2952 size_t counterIJ = m_CounterIJ[n];
2953
2954 // both species have a non-zero charge, and one is positive
2955 // and the other is negative
2956 if (charge(i)*charge(j) < 0.0) {
2957 m_BMX_IJ_LL[counterIJ] = m_Beta0MX_ij_LL[counterIJ]
2958 + m_Beta1MX_ij_LL[counterIJ] * m_gfunc_IJ[counterIJ]
2959 + m_Beta2MX_ij_LL[counterIJ] * m_g2func_IJ[counterIJ];
2960 if (Is > 1.0E-150) {
2961 m_BprimeMX_IJ_LL[counterIJ] = (m_Beta1MX_ij_LL[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
2962 m_Beta2MX_ij_LL[counterIJ] * m_h2func_IJ[counterIJ]/Is);
2963 } else {
2964 m_BprimeMX_IJ_LL[counterIJ] = 0.0;
2965 }
2966 m_BphiMX_IJ_LL[counterIJ] = m_BMX_IJ_LL[counterIJ] + Is*m_BprimeMX_IJ_LL[counterIJ];
2967 } else {
2968 m_BMX_IJ_LL[counterIJ] = 0.0;
2969 m_BprimeMX_IJ_LL[counterIJ] = 0.0;
2970 m_BphiMX_IJ_LL[counterIJ] = 0.0;
2971 }
2972 }
2973 }
2974
2975 // --------- SUBSECTION TO CALCULATE CMX_LL ----------
2976 for (size_t i = 1; i < m_kk-1; i++) {
2977 for (size_t j = i+1; j < m_kk; j++) {
2978 // Find the counterIJ for the symmetric binary interaction
2979 size_t n = m_kk*i + j;
2980 size_t counterIJ = m_CounterIJ[n];
2981
2982 // both species have a non-zero charge, and one is positive
2983 // and the other is negative
2984 if (charge(i)*charge(j) < 0.0) {
2985 m_CMX_IJ_LL[counterIJ] = m_CphiMX_ij_LL[counterIJ]/
2986 (2.0* sqrt(fabs(charge(i)*charge(j))));
2987 } else {
2988 m_CMX_IJ_LL[counterIJ] = 0.0;
2989 }
2990 }
2991 }
2992
2993 // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2994 for (size_t i = 1; i < m_kk-1; i++) {
2995 for (size_t j = i+1; j < m_kk; j++) {
2996 // Find the counterIJ for the symmetric binary interaction
2997 size_t n = m_kk*i + j;
2998 size_t counterIJ = m_CounterIJ[n];
2999
3000 // both species have a non-zero charge, and one is positive
3001 // and the other is negative
3002 if (charge(i)*charge(j) > 0) {
3003 m_Phi_IJ_LL[counterIJ] = m_Theta_ij_LL[counterIJ];
3004 m_Phiprime_IJ[counterIJ] = 0.0;
3005 m_PhiPhi_IJ_LL[counterIJ] = m_Phi_IJ_LL[counterIJ];
3006 } else {
3007 m_Phi_IJ_LL[counterIJ] = 0.0;
3008 m_Phiprime_IJ[counterIJ] = 0.0;
3009 m_PhiPhi_IJ_LL[counterIJ] = 0.0;
3010 }
3011 }
3012 }
3013
3014 // ----------- SUBSECTION FOR CALCULATION OF d2FdT2 ---------------------
3015 double d2AphidT2 = d2A_DebyedT2_TP() / 3.0;
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++) {
3020 // Find the counterIJ for the symmetric binary interaction
3021 size_t n = m_kk*i + j;
3022 size_t counterIJ = m_CounterIJ[n];
3023
3024 // both species have a non-zero charge, and one is positive
3025 // and the other is negative
3026 if (charge(i)*charge(j) < 0) {
3027 d2FdT2 += molality[i]*molality[j] * m_BprimeMX_IJ_LL[counterIJ];
3028 }
3029
3030 // Both species have a non-zero charge, and they
3031 // have the same sign, that is, both positive or both negative.
3032 if (charge(i)*charge(j) > 0) {
3033 d2FdT2 += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
3034 }
3035 }
3036 }
3037
3038 for (size_t i = 1; i < m_kk; i++) {
3039 // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
3040 if (charge(i) > 0) {
3041 // species i is the cation (positive) to calc the actcoeff
3042 double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3043 double sum1 = 0.0;
3044 double sum2 = 0.0;
3045 double sum3 = 0.0;
3046 double sum4 = 0.0;
3047 double sum5 = 0.0;
3048 for (size_t j = 1; j < m_kk; j++) {
3049 // Find the counterIJ for the symmetric binary interaction
3050 size_t n = m_kk*i + j;
3051 size_t counterIJ = m_CounterIJ[n];
3052
3053 if (charge(j) < 0.0) {
3054 // sum over all anions
3055 sum1 += molality[j]*
3056 (2.0*m_BMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3057 if (j < m_kk-1) {
3058 // This term is the ternary interaction involving the
3059 // non-duplicate sum over double anions, j, k, with
3060 // respect to the cation, i.
3061 for (size_t k = j+1; k < m_kk; k++) {
3062 // an inner sum over all anions
3063 if (charge(k) < 0.0) {
3064 n = k + j * m_kk + i * m_kk * m_kk;
3065 sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3066 }
3067 }
3068 }
3069 }
3070
3071 if (charge(j) > 0.0) {
3072 // sum over all cations
3073 if (j != i) {
3074 sum2 += molality[j]*(2.0*m_Phi_IJ_LL[counterIJ]);
3075 }
3076 for (size_t k = 1; k < m_kk; k++) {
3077 if (charge(k) < 0.0) {
3078 // two inner sums over anions
3079 n = k + j * m_kk + i * m_kk * m_kk;
3080 sum2 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3081
3082 // Find the counterIJ for the j,k interaction
3083 n = m_kk*j + k;
3084 size_t counterIJ2 = m_CounterIJ[n];
3085 sum4 += fabs(charge(i)) *
3086 molality[j]*molality[k]*m_CMX_IJ_LL[counterIJ2];
3087 }
3088 }
3089 }
3090
3091 // Handle neutral j species
3092 if (charge(j) == 0) {
3093 sum5 += molality[j]*2.0*m_Lambda_nj_LL(j,i);
3094 // Zeta interaction term
3095 for (size_t k = 1; k < m_kk; k++) {
3096 if (charge(k) < 0.0) {
3097 size_t izeta = j;
3098 size_t jzeta = i;
3099 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3100 double zeta_LL = m_Psi_ijk_LL[n];
3101 if (zeta_LL != 0.0) {
3102 sum5 += molality[j]*molality[k]*zeta_LL;
3103 }
3104 }
3105 }
3106 }
3107 }
3108 // Add all of the contributions up to yield the log of the
3109 // solute activity coefficients (molality scale)
3111 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3112 }
3113
3114 // ------ SUBSECTION FOR CALCULATING THE d2ACTCOEFFdT2 FOR ANIONS ------
3115 if (charge(i) < 0) {
3116 // species i is an anion (negative)
3117 double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3118 double sum1 = 0.0;
3119 double sum2 = 0.0;
3120 double sum3 = 0.0;
3121 double sum4 = 0.0;
3122 double sum5 = 0.0;
3123 for (size_t j = 1; j < m_kk; j++) {
3124 // Find the counterIJ for the symmetric binary interaction
3125 size_t n = m_kk*i + j;
3126 size_t counterIJ = m_CounterIJ[n];
3127
3128 // For Anions, do the cation interactions.
3129 if (charge(j) > 0) {
3130 sum1 += molality[j]*
3131 (2.0*m_BMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3132 if (j < m_kk-1) {
3133 for (size_t k = j+1; k < m_kk; k++) {
3134 // an inner sum over all cations
3135 if (charge(k) > 0) {
3136 n = k + j * m_kk + i * m_kk * m_kk;
3137 sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3138 }
3139 }
3140 }
3141 }
3142
3143 // For Anions, do the other anion interactions.
3144 if (charge(j) < 0.0) {
3145 // sum over all anions
3146 if (j != i) {
3147 sum2 += molality[j]*(2.0*m_Phi_IJ_LL[counterIJ]);
3148 }
3149 for (size_t k = 1; k < m_kk; k++) {
3150 if (charge(k) > 0.0) {
3151 // two inner sums over cations
3152 n = k + j * m_kk + i * m_kk * m_kk;
3153 sum2 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3154 // Find the counterIJ for the symmetric binary interaction
3155 n = m_kk*j + k;
3156 size_t counterIJ2 = m_CounterIJ[n];
3157 sum4 += fabs(charge(i)) *
3158 molality[j]*molality[k]*m_CMX_IJ_LL[counterIJ2];
3159 }
3160 }
3161 }
3162
3163 // for Anions, do the neutral species interaction
3164 if (charge(j) == 0.0) {
3165 sum5 += molality[j]*2.0*m_Lambda_nj_LL(j,i);
3166 // Zeta interaction term
3167 for (size_t k = 1; k < m_kk; k++) {
3168 if (charge(k) > 0.0) {
3169 size_t izeta = j;
3170 size_t jzeta = k;
3171 size_t kzeta = i;
3172 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3173 double zeta_LL = m_Psi_ijk_LL[n];
3174 if (zeta_LL != 0.0) {
3175 sum5 += molality[j]*molality[k]*zeta_LL;
3176 }
3177 }
3178 }
3179 }
3180 }
3182 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3183 }
3184
3185 // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
3186 // equations agree with my notes,
3187 // Equations agree with Pitzer,
3188 if (charge(i) == 0.0) {
3189 double sum1 = 0.0;
3190 double sum3 = 0.0;
3191 for (size_t j = 1; j < m_kk; j++) {
3192 sum1 += molality[j]*2.0*m_Lambda_nj_LL(i,j);
3193 // Zeta term -> we piggyback on the psi term
3194 if (charge(j) > 0.0) {
3195 for (size_t k = 1; k < m_kk; k++) {
3196 if (charge(k) < 0.0) {
3197 size_t n = k + j * m_kk + i * m_kk * m_kk;
3198 sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3199 }
3200 }
3201 }
3202 }
3203 double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i];
3204 m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3;
3205 }
3206 }
3207
3208 // ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 ---------
3209 double sum1 = 0.0;
3210 double sum2 = 0.0;
3211 double sum3 = 0.0;
3212 double sum4 = 0.0;
3213 double sum5 = 0.0;
3214 double sum6 = 0.0;
3215 double sum7 = 0.0;
3216
3217 // term1 is the temperature derivative of the DH term in the osmotic
3218 // coefficient expression
3219 // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
3220 // Is = Ionic strength on the molality scale (units of (gmol/kg))
3221 // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3222 double term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3223
3224 for (size_t j = 1; j < m_kk; j++) {
3225 // Loop Over Cations
3226 if (charge(j) > 0.0) {
3227 for (size_t k = 1; k < m_kk; k++) {
3228 if (charge(k) < 0.0) {
3229 // Find the counterIJ for the symmetric j,k binary interaction
3230 size_t n = m_kk*j + k;
3231 size_t counterIJ = m_CounterIJ[n];
3232
3233 sum1 += molality[j]*molality[k] *
3234 (m_BphiMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3235 }
3236 }
3237
3238 for (size_t k = j+1; k < m_kk; k++) {
3239 if (j == (m_kk-1)) {
3240 // we should never reach this step
3241 throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3242 "logic error 1 in Step 9 of hmw_act");
3243 }
3244 if (charge(k) > 0.0) {
3245 // Find the counterIJ for the symmetric j,k binary interaction
3246 // between 2 cations.
3247 size_t n = m_kk*j + k;
3248 size_t counterIJ = m_CounterIJ[n];
3249 sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_LL[counterIJ];
3250 for (size_t m = 1; m < m_kk; m++) {
3251 if (charge(m) < 0.0) {
3252 // species m is an anion
3253 n = m + k * m_kk + j * m_kk * m_kk;
3254 sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_LL[n];
3255 }
3256 }
3257 }
3258 }
3259 }
3260
3261 // Loop Over Anions
3262 if (charge(j) < 0) {
3263 for (size_t k = j+1; k < m_kk; k++) {
3264 if (j == m_kk-1) {
3265 // we should never reach this step
3266 throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3267 "logic error 2 in Step 9 of hmw_act");
3268 }
3269 if (charge(k) < 0) {
3270 // Find the counterIJ for the symmetric j,k binary interaction
3271 // between two anions
3272 size_t n = m_kk*j + k;
3273 size_t counterIJ = m_CounterIJ[n];
3274
3275 sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_LL[counterIJ];
3276 for (size_t m = 1; m < m_kk; m++) {
3277 if (charge(m) > 0.0) {
3278 n = m + k * m_kk + j * m_kk * m_kk;
3279 sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_LL[n];
3280 }
3281 }
3282 }
3283 }
3284 }
3285
3286 // Loop Over Neutral Species
3287 if (charge(j) == 0) {
3288 for (size_t k = 1; k < m_kk; k++) {
3289 if (charge(k) < 0.0) {
3290 sum4 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3291 }
3292 if (charge(k) > 0.0) {
3293 sum5 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3294 }
3295 if (charge(k) == 0.0) {
3296 if (k > j) {
3297 sum6 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3298 } else if (k == j) {
3299 sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3300 }
3301 }
3302 if (charge(k) < 0.0) {
3303 size_t izeta = j;
3304 for (size_t m = 1; m < m_kk; m++) {
3305 if (charge(m) > 0.0) {
3306 size_t jzeta = m;
3307 size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3308 double zeta_LL = m_Psi_ijk_LL[n];
3309 if (zeta_LL != 0.0) {
3310 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
3311 }
3312 }
3313 }
3314 }
3315 }
3316
3317 sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_LL[j];
3318 }
3319 }
3320 double sum_m_phi_minus_1 = 2.0 *
3321 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3322 // Calculate the osmotic coefficient from
3323 // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
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);
3327 } else {
3328 d2_osmotic_coef_dT2 = 0.0;
3329 }
3330 double d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
3331
3332 // In Cantera, we define the activity coefficient of the solvent as
3333 //
3334 // act_0 = actcoeff_0 * Xmol_0
3335 //
3336 // We have just computed act_0. However, this routine returns
3337 // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3338 m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2;
3339}
3340
3342{
3343 static const int cacheId = m_cache.getId();
3344 CachedScalar cached = m_cache.getScalar(cacheId);
3345 if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
3346 return;
3347 }
3348
3351
3352 for (size_t k = 1; k < m_kk; k++) {
3353 if (CROP_speciesCropped_[k] == 2) {
3355 }
3356 }
3357
3358 if (CROP_speciesCropped_[0]) {
3360 }
3361
3363}
3364
3366{
3367 const double* molality = m_molalitiesCropped.data();
3368
3369 // Local variables defined by Coltrin
3370 double etheta[5][5], etheta_prime[5][5], sqrtIs;
3371
3372 // Molality based ionic strength of the solution
3373 double Is = 0.0;
3374
3375 // Molar charge of the solution: In Pitzer's notation, this is his variable
3376 // called "Z".
3377 double molarcharge = 0.0;
3378
3379 // molalitysum is the sum of the molalities over all solutes, even those
3380 // with zero charge.
3381 double molalitysum = 0.0;
3382 double currTemp = temperature();
3383 double currPres = pressure();
3384
3385 // Make sure the counter variables are setup
3387
3388 // ---------- Calculate common sums over solutes ---------------------
3389 for (size_t n = 1; n < m_kk; n++) {
3390 // ionic strength
3391 Is += charge(n) * charge(n) * molality[n];
3392 // total molar charge
3393 molarcharge += fabs(charge(n)) * molality[n];
3394 molalitysum += molality[n];
3395 }
3396 Is *= 0.5;
3397
3398 // Store the ionic molality in the object for reference.
3399 m_IionicMolality = Is;
3400 sqrtIs = sqrt(Is);
3401
3402 // The following call to calc_lambdas() calculates all 16 elements of the
3403 // elambda and elambda1 arrays, given the value of the ionic strength (Is)
3404 calc_lambdas(Is);
3405
3406 // Step 2: Find the coefficients E-theta and E-thetaprime for all
3407 // combinations of positive unlike charges up to 4
3408 for (int z1 = 1; z1 <=4; z1++) {
3409 for (int z2 =1; z2 <=4; z2++) {
3410 calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
3411 }
3412 }
3413
3414 // calculate g(x) and hfunc(x) for each cation-anion pair MX
3415 // In the original literature, hfunc, was called gprime. However,
3416 // it's not the derivative of g(x), so I renamed it.
3417 for (size_t i = 1; i < (m_kk - 1); i++) {
3418 for (size_t j = (i+1); j < m_kk; j++) {
3419 // Find the counterIJ for the symmetric binary interaction
3420 size_t n = m_kk*i + j;
3421 size_t counterIJ = m_CounterIJ[n];
3422
3423 // Only loop over oppositely charge species
3424 if (charge(i)*charge(j) < 0) {
3425 // x is a reduced function variable
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);
3429 m_hfunc_IJ[counterIJ] = -2.0*
3430 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
3431 } else {
3432 m_gfunc_IJ[counterIJ] = 0.0;
3433 m_hfunc_IJ[counterIJ] = 0.0;
3434 }
3435
3436 if (m_Beta2MX_ij_P[counterIJ] != 0.0) {
3437 double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
3438 if (x2 > 1.0E-100) {
3439 m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3440 m_h2func_IJ[counterIJ] = -2.0 *
3441 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3442 } else {
3443 m_g2func_IJ[counterIJ] = 0.0;
3444 m_h2func_IJ[counterIJ] = 0.0;
3445 }
3446 }
3447 } else {
3448 m_gfunc_IJ[counterIJ] = 0.0;
3449 m_hfunc_IJ[counterIJ] = 0.0;
3450 }
3451 }
3452 }
3453
3454 // SUBSECTION TO CALCULATE BMX_P, BprimeMX_P, BphiMX_P
3455 // These are now temperature derivatives of the previously calculated
3456 // quantities.
3457 for (size_t i = 1; i < m_kk - 1; i++) {
3458 for (size_t j = i+1; j < m_kk; j++) {
3459 // Find the counterIJ for the symmetric binary interaction
3460 size_t n = m_kk*i + j;
3461 size_t counterIJ = m_CounterIJ[n];
3462
3463 // both species have a non-zero charge, and one is positive
3464 // and the other is negative
3465 if (charge(i)*charge(j) < 0.0) {
3466 m_BMX_IJ_P[counterIJ] = m_Beta0MX_ij_P[counterIJ]
3467 + m_Beta1MX_ij_P[counterIJ] * m_gfunc_IJ[counterIJ]
3468 + m_Beta2MX_ij_P[counterIJ] * m_g2func_IJ[counterIJ];
3469 if (Is > 1.0E-150) {
3470 m_BprimeMX_IJ_P[counterIJ] = (m_Beta1MX_ij_P[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
3471 m_Beta2MX_ij_P[counterIJ] * m_h2func_IJ[counterIJ]/Is);
3472 } else {
3473 m_BprimeMX_IJ_P[counterIJ] = 0.0;
3474 }
3475 m_BphiMX_IJ_P[counterIJ] = m_BMX_IJ_P[counterIJ] + Is*m_BprimeMX_IJ_P[counterIJ];
3476 } else {
3477 m_BMX_IJ_P[counterIJ] = 0.0;
3478 m_BprimeMX_IJ_P[counterIJ] = 0.0;
3479 m_BphiMX_IJ_P[counterIJ] = 0.0;
3480 }
3481 }
3482 }
3483
3484 // --------- SUBSECTION TO CALCULATE CMX_P ----------
3485 for (size_t i = 1; i < m_kk-1; i++) {
3486 for (size_t j = i+1; j < m_kk; j++) {
3487 // Find the counterIJ for the symmetric binary interaction
3488 size_t n = m_kk*i + j;
3489 size_t counterIJ = m_CounterIJ[n];
3490
3491 // both species have a non-zero charge, and one is positive
3492 // and the other is negative
3493 if (charge(i)*charge(j) < 0.0) {
3494 m_CMX_IJ_P[counterIJ] = m_CphiMX_ij_P[counterIJ]/
3495 (2.0* sqrt(fabs(charge(i)*charge(j))));
3496 } else {
3497 m_CMX_IJ_P[counterIJ] = 0.0;
3498 }
3499 }
3500 }
3501
3502 // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
3503 for (size_t i = 1; i < m_kk-1; i++) {
3504 for (size_t j = i+1; j < m_kk; j++) {
3505 // Find the counterIJ for the symmetric binary interaction
3506 size_t n = m_kk*i + j;
3507 size_t counterIJ = m_CounterIJ[n];
3508
3509 // both species have a non-zero charge, and one is positive
3510 // and the other is negative
3511 if (charge(i)*charge(j) > 0) {
3512 m_Phi_IJ_P[counterIJ] = m_Theta_ij_P[counterIJ];
3513 m_Phiprime_IJ[counterIJ] = 0.0;
3514 m_PhiPhi_IJ_P[counterIJ] = m_Phi_IJ_P[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
3515 } else {
3516 m_Phi_IJ_P[counterIJ] = 0.0;
3517 m_Phiprime_IJ[counterIJ] = 0.0;
3518 m_PhiPhi_IJ_P[counterIJ] = 0.0;
3519 }
3520 }
3521 }
3522
3523 // ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
3524 double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres);
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++) {
3530 // Find the counterIJ for the symmetric binary interaction
3531 size_t n = m_kk*i + j;
3532 size_t counterIJ = m_CounterIJ[n];
3533
3534 // both species have a non-zero charge, and one is positive
3535 // and the other is negative
3536 if (charge(i)*charge(j) < 0) {
3537 dFdP += molality[i]*molality[j] * m_BprimeMX_IJ_P[counterIJ];
3538 }
3539
3540 // Both species have a non-zero charge, and they
3541 // have the same sign, that is, both positive or both negative.
3542 if (charge(i)*charge(j) > 0) {
3543 dFdP += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
3544 }
3545 }
3546 }
3547
3548 for (size_t i = 1; i < m_kk; i++) {
3549 // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS -----
3550 if (charge(i) > 0) {
3551 // species i is the cation (positive) to calc the actcoeff
3552 double zsqdFdP = charge(i)*charge(i)*dFdP;
3553 double sum1 = 0.0;
3554 double sum2 = 0.0;
3555 double sum3 = 0.0;
3556 double sum4 = 0.0;
3557 double sum5 = 0.0;
3558 for (size_t j = 1; j < m_kk; j++) {
3559 // Find the counterIJ for the symmetric binary interaction
3560 size_t n = m_kk*i + j;
3561 size_t counterIJ = m_CounterIJ[n];
3562
3563 if (charge(j) < 0.0) {
3564 // sum over all anions
3565 sum1 += molality[j]*
3566 (2.0*m_BMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3567 if (j < m_kk-1) {
3568 // This term is the ternary interaction involving the
3569 // non-duplicate sum over double anions, j, k, with
3570 // respect to the cation, i.
3571 for (size_t k = j+1; k < m_kk; k++) {
3572 // an inner sum over all anions
3573 if (charge(k) < 0.0) {
3574 n = k + j * m_kk + i * m_kk * m_kk;
3575 sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3576 }
3577 }
3578 }
3579 }
3580
3581 if (charge(j) > 0.0) {
3582 // sum over all cations
3583 if (j != i) {
3584 sum2 += molality[j]*(2.0*m_Phi_IJ_P[counterIJ]);
3585 }
3586 for (size_t k = 1; k < m_kk; k++) {
3587 if (charge(k) < 0.0) {
3588 // two inner sums over anions
3589 n = k + j * m_kk + i * m_kk * m_kk;
3590 sum2 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3591
3592 // Find the counterIJ for the j,k interaction
3593 n = m_kk*j + k;
3594 size_t counterIJ2 = m_CounterIJ[n];
3595 sum4 += fabs(charge(i)) *
3596 molality[j]*molality[k]*m_CMX_IJ_P[counterIJ2];
3597 }
3598 }
3599 }
3600
3601 // for Anions, do the neutral species interaction
3602 if (charge(j) == 0) {
3603 sum5 += molality[j]*2.0*m_Lambda_nj_P(j,i);
3604 // Zeta interaction term
3605 for (size_t k = 1; k < m_kk; k++) {
3606 if (charge(k) < 0.0) {
3607 size_t izeta = j;
3608 size_t jzeta = i;
3609 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3610 double zeta_P = m_Psi_ijk_P[n];
3611 if (zeta_P != 0.0) {
3612 sum5 += molality[j]*molality[k]*zeta_P;
3613 }
3614 }
3615 }
3616 }
3617 }
3618
3619 // Add all of the contributions up to yield the log of the
3620 // solute activity coefficients (molality scale)
3622 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3623 }
3624
3625 // ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR ANIONS ------
3626 if (charge(i) < 0) {
3627 // species i is an anion (negative)
3628 double zsqdFdP = charge(i)*charge(i)*dFdP;
3629 double sum1 = 0.0;
3630 double sum2 = 0.0;
3631 double sum3 = 0.0;
3632 double sum4 = 0.0;
3633 double sum5 = 0.0;
3634 for (size_t j = 1; j < m_kk; j++) {
3635 // Find the counterIJ for the symmetric binary interaction
3636 size_t n = m_kk*i + j;
3637 size_t counterIJ = m_CounterIJ[n];
3638
3639 // For Anions, do the cation interactions.
3640 if (charge(j) > 0) {
3641 sum1 += molality[j] *
3642 (2.0*m_BMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3643 if (j < m_kk-1) {
3644 for (size_t k = j+1; k < m_kk; k++) {
3645 // an inner sum over all cations
3646 if (charge(k) > 0) {
3647 n = k + j * m_kk + i * m_kk * m_kk;
3648 sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3649 }
3650 }
3651 }
3652 }
3653
3654 // For Anions, do the other anion interactions.
3655 if (charge(j) < 0.0) {
3656 // sum over all anions
3657 if (j != i) {
3658 sum2 += molality[j]*(2.0*m_Phi_IJ_P[counterIJ]);
3659 }
3660 for (size_t k = 1; k < m_kk; k++) {
3661 if (charge(k) > 0.0) {
3662 // two inner sums over cations
3663 n = k + j * m_kk + i * m_kk * m_kk;
3664 sum2 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3665 // Find the counterIJ for the symmetric binary interaction
3666 n = m_kk*j + k;
3667 size_t counterIJ2 = m_CounterIJ[n];
3668 sum4 += fabs(charge(i))*
3669 molality[j]*molality[k]*m_CMX_IJ_P[counterIJ2];
3670 }
3671 }
3672 }
3673
3674 // for Anions, do the neutral species interaction
3675 if (charge(j) == 0.0) {
3676 sum5 += molality[j]*2.0*m_Lambda_nj_P(j,i);
3677 // Zeta interaction term
3678 for (size_t k = 1; k < m_kk; k++) {
3679 if (charge(k) > 0.0) {
3680 size_t izeta = j;
3681 size_t jzeta = k;
3682 size_t kzeta = i;
3683 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3684 double zeta_P = m_Psi_ijk_P[n];
3685 if (zeta_P != 0.0) {
3686 sum5 += molality[j]*molality[k]*zeta_P;
3687 }
3688 }
3689 }
3690 }
3691 }
3693 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3694 }
3695
3696 // ------ SUBSECTION FOR CALCULATING d NEUTRAL SOLUTE ACT COEFF dP -----
3697 if (charge(i) == 0.0) {
3698 double sum1 = 0.0;
3699 double sum3 = 0.0;
3700 for (size_t j = 1; j < m_kk; j++) {
3701 sum1 += molality[j]*2.0*m_Lambda_nj_P(i,j);
3702 // Zeta term -> we piggyback on the psi term
3703 if (charge(j) > 0.0) {
3704 for (size_t k = 1; k < m_kk; k++) {
3705 if (charge(k) < 0.0) {
3706 size_t n = k + j * m_kk + i * m_kk * m_kk;
3707 sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3708 }
3709 }
3710 }
3711 }
3712 double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i];
3713 m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3;
3714 }
3715 }
3716
3717 // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP ---------
3718 double sum1 = 0.0;
3719 double sum2 = 0.0;
3720 double sum3 = 0.0;
3721 double sum4 = 0.0;
3722 double sum5 = 0.0;
3723 double sum6 = 0.0;
3724 double sum7 = 0.0;
3725
3726 // term1 is the temperature derivative of the DH term in the osmotic
3727 // coefficient expression
3728 // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
3729 // Is = Ionic strength on the molality scale (units of (gmol/kg))
3730 // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3731 double term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3732
3733 for (size_t j = 1; j < m_kk; j++) {
3734 // Loop Over Cations
3735 if (charge(j) > 0.0) {
3736 for (size_t k = 1; k < m_kk; k++) {
3737 if (charge(k) < 0.0) {
3738 // Find the counterIJ for the symmetric j,k binary interaction
3739 size_t n = m_kk*j + k;
3740 size_t counterIJ = m_CounterIJ[n];
3741 sum1 += molality[j]*molality[k]*
3742 (m_BphiMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3743 }
3744 }
3745
3746 for (size_t k = j+1; k < m_kk; k++) {
3747 if (j == (m_kk-1)) {
3748 // we should never reach this step
3749 throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3750 "logic error 1 in Step 9 of hmw_act");
3751 }
3752 if (charge(k) > 0.0) {
3753 // Find the counterIJ for the symmetric j,k binary interaction
3754 // between 2 cations.
3755 size_t n = m_kk*j + k;
3756 size_t counterIJ = m_CounterIJ[n];
3757 sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_P[counterIJ];
3758 for (size_t m = 1; m < m_kk; m++) {
3759 if (charge(m) < 0.0) {
3760 // species m is an anion
3761 n = m + k * m_kk + j * m_kk * m_kk;
3762 sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_P[n];
3763 }
3764 }
3765 }
3766 }
3767 }
3768
3769 // Loop Over Anions
3770 if (charge(j) < 0) {
3771 for (size_t k = j+1; k < m_kk; k++) {
3772 if (j == m_kk-1) {
3773 // we should never reach this step
3774 throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3775 "logic error 2 in Step 9 of hmw_act");
3776 }
3777 if (charge(k) < 0) {
3778 // Find the counterIJ for the symmetric j,k binary interaction
3779 // between two anions
3780 size_t n = m_kk*j + k;
3781 size_t counterIJ = m_CounterIJ[n];
3782
3783 sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_P[counterIJ];
3784 for (size_t m = 1; m < m_kk; m++) {
3785 if (charge(m) > 0.0) {
3786 n = m + k * m_kk + j * m_kk * m_kk;
3787 sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_P[n];
3788 }
3789 }
3790 }
3791 }
3792 }
3793
3794 // Loop Over Neutral Species
3795 if (charge(j) == 0) {
3796 for (size_t k = 1; k < m_kk; k++) {
3797 if (charge(k) < 0.0) {
3798 sum4 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3799 }
3800 if (charge(k) > 0.0) {
3801 sum5 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3802 }
3803 if (charge(k) == 0.0) {
3804 if (k > j) {
3805 sum6 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3806 } else if (k == j) {
3807 sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3808 }
3809 }
3810 if (charge(k) < 0.0) {
3811 size_t izeta = j;
3812 for (size_t m = 1; m < m_kk; m++) {
3813 if (charge(m) > 0.0) {
3814 size_t jzeta = m;
3815 size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3816 double zeta_P = m_Psi_ijk_P[n];
3817 if (zeta_P != 0.0) {
3818 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
3819 }
3820 }
3821 }
3822 }
3823 }
3824
3825 sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_P[j];
3826 }
3827 }
3828 double sum_m_phi_minus_1 = 2.0 *
3829 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3830
3831 // Calculate the osmotic coefficient from
3832 // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
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);
3836 } else {
3837 d_osmotic_coef_dP = 0.0;
3838 }
3839 double d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
3840
3841 // In Cantera, we define the activity coefficient of the solvent as
3842 //
3843 // act_0 = actcoeff_0 * Xmol_0
3844 //
3845 // We have just computed act_0. However, this routine returns
3846 // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3847 m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP;
3848}
3849
3850void HMWSoln::calc_lambdas(double is) const
3851{
3852 if( m_last_is == is ) {
3853 return;
3854 }
3855 m_last_is = is;
3856
3857 // Coefficients c1-c4 are used to approximate the integral function "J";
3858 // aphi is the Debye-Huckel constant at 25 C
3859 double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
3860 double aphi = 0.392; /* Value at 25 C */
3861 if (is < 1.0E-150) {
3862 for (int i = 0; i < 17; i++) {
3863 elambda[i] = 0.0;
3864 elambda1[i] = 0.0;
3865 }
3866 return;
3867 }
3868
3869 // Calculate E-lambda terms for charge combinations of like sign,
3870 // using method of Pitzer (1975). Charges up to 4 are calculated.
3871 for (int i=1; i<=4; i++) {
3872 for (int j=i; j<=4; j++) {
3873 int ij = i*j;
3874
3875 // calculate the product of the charges
3876 double zprod = (double)ij;
3877
3878 // calculate Xmn (A1) from Harvie, Weare (1980).
3879 double x = 6.0* zprod * aphi * sqrt(is); // eqn 23
3880
3881 double jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4))); // eqn 47
3882
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);
3886
3887 elambda[ij] = zprod*jfunc / (4.0*is); // eqn 14
3888 elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
3889 - elambda[ij])/is;
3890 }
3891 }
3892}
3893
3894void HMWSoln::calc_thetas(int z1, int z2,
3895 double* etheta, double* etheta_prime) const
3896{
3897 // Calculate E-theta(i) and E-theta'(I) using method of Pitzer (1987)
3898 int i = abs(z1);
3899 int j = abs(z2);
3900
3901 AssertThrowMsg(i <= 4 && j <= 4, "HMWSoln::calc_thetas",
3902 "we shouldn't be here");
3903 AssertThrowMsg(i != 0 && j != 0, "HMWSoln::calc_thetas",
3904 "called with one species being neutral");
3905
3906 // Check to see if the charges are of opposite sign. If they are of opposite
3907 // sign then their etheta interaction is zero.
3908 if (z1*z2 < 0) {
3909 *etheta = 0.0;
3910 *etheta_prime = 0.0;
3911 } else {
3912 // Actually calculate the interaction.
3913 double f1 = (double)i / (2.0 * j);
3914 double f2 = (double)j / (2.0 * i);
3915 *etheta = elambda[i*j] - f1*elambda[j*j] - f2*elambda[i*i];
3916 *etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i];
3917 }
3918}
3919
3921{
3922 // Calculate the molalities. Currently, the molalities may not be current
3923 // with respect to the contents of the State objects' data.
3925 double xmolSolvent = moleFraction(0);
3926 double xx = std::max(m_xmolSolventMIN, xmolSolvent);
3927 // Exponentials - trial 2
3928 if (xmolSolvent > IMS_X_o_cutoff_) {
3929 for (size_t k = 1; k < m_kk; k++) {
3930 IMS_lnActCoeffMolal_[k]= 0.0;
3931 }
3932 IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
3933 return;
3934 } else {
3935 double xoverc = xmolSolvent/IMS_cCut_;
3936 double eterm = std::exp(-xoverc);
3937
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));
3943
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));
3949
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);
3953
3954 tmp = log(xx) + lngammak;
3955 for (size_t k = 1; k < m_kk; k++) {
3956 IMS_lnActCoeffMolal_[k]= tmp;
3957 }
3958 IMS_lnActCoeffMolal_[0] = lngammao;
3959 }
3960 return;
3961}
3962
3964{
3966 vector<double>& moleF = m_tmpV;
3967
3968 // Update the coefficients wrt Temperature. Calculate the derivatives as well
3970 getMoleFractions(moleF.data());
3971
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",
3975 k, speciesName(k), moleF[k], m_molalitiesCropped[k], charge(k));
3976 }
3977
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;
3983 size_t ct = m_CounterIJ[n];
3984 writelogf(" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
3985 speciesName(i), speciesName(j),
3986 m_Beta0MX_ij[ct], m_Beta1MX_ij[ct],
3987 m_Beta2MX_ij[ct], m_CphiMX_ij[ct],
3988 m_Alpha1MX_ij[ct], m_Theta_ij[ct]);
3989 }
3990 }
3991
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++) {
3996 size_t n = k + j * m_kk + i * m_kk * m_kk;
3997 if (m_Psi_ijk[n] != 0.0) {
3998 writelogf(" %-16s %-16s %-16s %9.5f \n",
3999 speciesName(i), speciesName(j),
4000 speciesName(k), m_Psi_ijk[n]);
4001 }
4002 }
4003 }
4004 }
4005}
4006
4007void HMWSoln::applyphScale(double* acMolality) const
4008{
4010 return;
4011 }
4013 double lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
4014 double lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
4015 double afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4016 for (size_t k = 0; k < m_kk; k++) {
4017 acMolality[k] *= exp(charge(k) * afac);
4018 }
4019}
4020
4022{
4025 return;
4026 }
4028 double lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
4029 double lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
4030 double afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4031 for (size_t k = 0; k < m_kk; k++) {
4033 }
4034}
4035
4037{
4040 return;
4041 }
4043 double dlnGammaClM_dT_s2 = s_NBS_CLM_dlnMolalityActCoeff_dT();
4044 double dlnGammaCLM_dT_s1 = m_dlnActCoeffMolaldT_Unscaled[m_indexCLM];
4045 double afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
4046 for (size_t k = 0; k < m_kk; k++) {
4048 }
4049}
4050
4052{
4055 return;
4056 }
4058 double d2lnGammaClM_dT2_s2 = s_NBS_CLM_d2lnMolalityActCoeff_dT2();
4059 double d2lnGammaCLM_dT2_s1 = m_d2lnActCoeffMolaldT2_Unscaled[m_indexCLM];
4060 double afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
4061 for (size_t k = 0; k < m_kk; k++) {
4063 }
4064}
4065
4067{
4070 return;
4071 }
4073 double dlnGammaClM_dP_s2 = s_NBS_CLM_dlnMolalityActCoeff_dP();
4074 double dlnGammaCLM_dP_s1 = m_dlnActCoeffMolaldP_Unscaled[m_indexCLM];
4075 double afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
4076 for (size_t k = 0; k < m_kk; k++) {
4078 }
4079}
4080
4082{
4083 double sqrtIs = sqrt(m_IionicMolality);
4084 double A = A_Debye_TP();
4085 double lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
4086 return lnGammaClMs2;
4087}
4088
4090{
4091 double sqrtIs = sqrt(m_IionicMolality);
4092 double dAdT = dA_DebyedT_TP();
4093 return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
4094}
4095
4097{
4098 double sqrtIs = sqrt(m_IionicMolality);
4099 double d2AdT2 = d2A_DebyedT2_TP();
4100 return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
4101}
4102
4104{
4105 double sqrtIs = sqrt(m_IionicMolality);
4106 double dAdP = dA_DebyedP_TP();
4107 return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
4108}
4109
4110}
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.
Definition AnyMap.h:427
size_t size() const
Returns the number of elements in this map.
Definition AnyMap.h:622
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition Array.h:203
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition Array.cpp:47
Base class for exceptions thrown by Cantera classes.
void calcIMSCutoffParams_()
Precalculate the IMS Cutoff parameters for typeCutoff = 2.
Definition HMWSoln.cpp:1473
vector< double > m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT.
Definition HMWSoln.h:1721
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition HMWSoln.h:1575
vector< double > m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
Definition HMWSoln.h:1668
vector< double > m_Beta1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition HMWSoln.h:1459
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
Definition HMWSoln.cpp:54
void applyphScale(double *acMolality) const override
Apply the current phScale to a set of activity Coefficients or activities.
Definition HMWSoln.cpp:4007
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition HMWSoln.h:1423
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer's papers.
Definition HMWSoln.h:1452
double m_A_Debye
A_Debye: this expression appears on the top of the ln actCoeff term in the general Debye-Huckel expre...
Definition HMWSoln.h:1414
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...
Definition HMWSoln.cpp:1104
vector< double > m_h2func_IJ
hfunc2, was called gprime in Pitzer's paper.
Definition HMWSoln.h:1772
vector< double > m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1821
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
Definition HMWSoln.cpp:231
vector< double > m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition HMWSoln.h:1435
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
Definition HMWSoln.cpp:211
double IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
Definition HMWSoln.h:1876
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
Definition HMWSoln.h:1635
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature.
Definition HMWSoln.cpp:4036
vector< double > m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
Definition HMWSoln.h:1761
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
Definition HMWSoln.h:1358
void calcMCCutoffParams_()
Calculate molality cut-off parameters.
Definition HMWSoln.cpp:1526
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer's papers.
Definition HMWSoln.h:1502
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition HMWSoln.cpp:2830
vector< double > m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1492
vector< double > m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1441
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
Definition HMWSoln.h:1382
vector< double > m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
Definition HMWSoln.h:1735
HMWSoln(const string &inputFile="", const string &id="")
Construct and initialize an HMWSoln ThermoPhase object directly from an input file.
Definition HMWSoln.cpp:41
vector< double > m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
Definition HMWSoln.h:1698
vector< double > m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1533
vector< double > m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1835
vector< double > m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1560
vector< double > m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1557
vector< double > m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1530
vector< double > m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1838
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
Definition HMWSoln.h:1629
double s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale.
Definition HMWSoln.cpp:4089
vector< double > m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
Definition HMWSoln.h:1592
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
Definition HMWSoln.cpp:752
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
Definition HMWSoln.cpp:576
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure.
Definition HMWSoln.cpp:4066
vector< double > m_tmpV
vector of size m_kk, used as a temporary holding area.
Definition HMWSoln.h:1426
void printCoeffs() const
Print out all of the input Pitzer coefficients.
Definition HMWSoln.cpp:3963
void getActivityConcentrations(double *c) const override
This method returns an array of generalized activity concentrations.
Definition HMWSoln.cpp:159
vector< double > m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
Definition HMWSoln.h:1588
vector< double > m_hfunc_IJ
hfunc, was called gprime in Pitzer's paper.
Definition HMWSoln.h:1768
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
Definition HMWSoln.cpp:2311
void s_update_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition HMWSoln.cpp:1245
vector< double > m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1468
vector< int > CROP_speciesCropped_
This is a boolean-type vector indicating whether a species's activity coefficient is in the cropped r...
Definition HMWSoln.h:1915
vector< double > m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition HMWSoln.h:1517
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_V.
Definition HMWSoln.cpp:1081
vector< double > m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1805
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
Definition HMWSoln.cpp:285
vector< double > m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1824
vector< double > m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1811
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
Definition HMWSoln.cpp:1450
vector< double > m_BMX_IJ
Intermediate variable called BMX in Pitzer's paper.
Definition HMWSoln.h:1776
vector< double > m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P.
Definition HMWSoln.h:1729
vector< double > m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1851
vector< double > m_BMX_IJ_P
Derivative of BMX_IJ wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1785
double cv_mole() const override
Molar heat capacity at constant volume. Units: J/kmol/K.
Definition HMWSoln.cpp:130
vector< double > m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
Definition HMWSoln.h:1688
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
Definition HMWSoln.h:1616
void getUnscaledMolalityActivityCoefficients(double *acMolality) const override
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
Definition HMWSoln.cpp:198
void setA_Debye(double A)
Set the A_Debye parameter.
Definition HMWSoln.cpp:539
vector< double > IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
Definition HMWSoln.h:1862
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model.
Definition HMWSoln.h:1353
vector< double > m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T.
Definition HMWSoln.h:1713
vector< double > m_BMX_IJ_L
Derivative of BMX_IJ wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1779
virtual double relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
Definition HMWSoln.cpp:60
vector< double > m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1438
vector< double > m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1444
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients' dependence on the temperature.
Definition HMWSoln.cpp:1562
double s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
Definition HMWSoln.cpp:4081
vector< double > m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition HMWSoln.h:1554
vector< double > m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1486
vector< double > m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
Definition HMWSoln.h:1678
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients.
Definition HMWSoln.cpp:3365
vector< double > m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer's paper.
Definition HMWSoln.h:1802
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
Definition HMWSoln.h:1626
vector< double > m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer's paper.
Definition HMWSoln.h:1789
vector< double > m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1563
vector< double > m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer's paper.
Definition HMWSoln.h:1832
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
Definition HMWSoln.cpp:3894
vector< int > m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species.
Definition HMWSoln.h:1746
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
Definition HMWSoln.h:1691
virtual double relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
Definition HMWSoln.cpp:72
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_L.
Definition HMWSoln.cpp:1070
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
Definition HMWSoln.cpp:112
PDSS * m_waterSS
Water standard state calculator.
Definition HMWSoln.h:1420
void calcDensity() override
Calculate the density of the mixture using the partial molar volumes and mole fractions as input.
Definition HMWSoln.cpp:142
void calc_lambdas(double is) const
Calculate the lambda interactions.
Definition HMWSoln.cpp:3850
double elambda1[17]
This is elambda1, MEC.
Definition HMWSoln.h:1752
vector< double > m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1465
vector< double > m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
Definition HMWSoln.h:1705
double s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
Definition HMWSoln.cpp:4103
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition HMWSoln.h:1362
vector< double > m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T.
Definition HMWSoln.h:1709
vector< double > m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition HMWSoln.h:1483
vector< double > m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1795
vector< double > m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer's paper.
Definition HMWSoln.h:1828
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients.
Definition HMWSoln.cpp:2339
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
Definition HMWSoln.cpp:1808
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
Definition HMWSoln.cpp:124
vector< double > m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1841
void getActivities(double *ac) const override
Get the array of non-dimensional activities at the current solution temperature, pressure,...
Definition HMWSoln.cpp:182
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...
Definition HMWSoln.cpp:1014
vector< double > m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1792
vector< double > m_Phi_IJ
Intermediate variable called Phi in Pitzer's paper.
Definition HMWSoln.h:1815
vector< double > m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1782
vector< double > m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT.
Definition HMWSoln.h:1717
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
Definition HMWSoln.h:1632
double IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process;
Definition HMWSoln.h:1866
double IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
Definition HMWSoln.h:1869
vector< double > m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1808
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition HMWSoln.cpp:2858
vector< double > m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1848
void getPartialMolarCp(double *cpbar) const override
Return an array of partial molar heat capacities for the species in the mixture.
Definition HMWSoln.cpp:298
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
Definition HMWSoln.h:1650
vector< double > m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
Definition HMWSoln.h:1596
void initLengths()
Initialize all of the species-dependent lengths in the object.
Definition HMWSoln.cpp:1132
vector< double > m_g2func_IJ
This is the value of g2(x2) in Pitzer's papers. Vector index is counterIJ.
Definition HMWSoln.h:1764
double gibbs_mole() const override
Molar Gibbs function. Units: J/kmol.
Definition HMWSoln.cpp:118
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
Definition HMWSoln.cpp:172
vector< double > m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition HMWSoln.h:1524
vector< double > m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1798
vector< double > m_CMX_IJ
Intermediate variable called CMX in Pitzer's paper.
Definition HMWSoln.h:1845
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
Definition HMWSoln.cpp:3341
double MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process;
Definition HMWSoln.h:1892
double satPressure(double T) override
Get the saturation pressure for a given temperature.
Definition HMWSoln.cpp:318
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
Definition HMWSoln.h:1365
vector< double > m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT. Vector index is counterIJ.
Definition HMWSoln.h:1489
double elambda[17]
This is elambda, MEC.
Definition HMWSoln.h:1749
vector< double > m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P.
Definition HMWSoln.h:1725
vector< double > m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
Definition HMWSoln.h:1658
vector< double > m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
Definition HMWSoln.h:1600
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition HMWSoln.cpp:3920
void calcMolalitiesCropped() const
Calculate the cropped molalities.
Definition HMWSoln.cpp:1309
vector< double > m_CphiMX_ij_L
Derivative of Cphi_ij[i][j] wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1527
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_J.
Definition HMWSoln.cpp:1092
void s_updateScaling_pHScaling_dT2() const
Apply the current phScale to a set of 2nd derivatives of the activity Coefficients wrt temperature.
Definition HMWSoln.cpp:4051
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
Definition HMWSoln.cpp:250
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer's papers.
Definition HMWSoln.h:1476
vector< double > m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P. Vector index is counterIJ.
Definition HMWSoln.h:1854
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
Definition HMWSoln.cpp:4021
vector< double > m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1462
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
Definition HMWSoln.h:1542
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.
Definition HMWSoln.cpp:982
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...
Definition HMWSoln.cpp:1038
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
Definition HMWSoln.h:1738
vector< double > m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T. Vector index is counterIJ.
Definition HMWSoln.h:1818
vector< double > m_gamma_tmp
Intermediate storage of the activity coefficient itself.
Definition HMWSoln.h:1858
double s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale.
Definition HMWSoln.cpp:4096
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition AnyMap.h:738
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.
Definition PDSS_Water.h:50
virtual double satPressure(double T)
saturation pressure
Definition PDSS.cpp:157
virtual void setState_TP(double temp, double pres)
Set the internal temperature and pressure.
Definition PDSS.cpp:152
void assignDensity(const double density_)
Set the internally stored constant density (kg/m^3) of the phase.
Definition Phase.cpp:597
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
Definition Phase.cpp:289
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
Definition Phase.h:822
size_t m_kk
Number of species in the phase.
Definition Phase.h:842
double temperature() const
Temperature (K).
Definition Phase.h:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:655
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:434
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:129
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition Phase.cpp:439
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition Phase.h:761
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition Phase.cpp:616
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition Phase.cpp:856
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition Phase.cpp:581
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition Phase.h:538
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.
Definition ValueCache.h:161
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.
Definition global.h:195
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
int sign(T x)
Sign of a number.
Definition global.h:336
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
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.
Definition ct_defs.h:158
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.
Definition ValueCache.h:33
T value
The value of the cached property.
Definition ValueCache.h:113
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.
Definition ValueCache.h:39