Cantera  3.1.0a1
Loading...
Searching...
No Matches
DebyeHuckel.cpp
Go to the documentation of this file.
1/**
2 * @file DebyeHuckel.cpp
3 * Declarations for the DebyeHuckel ThermoPhase object, which models dilute
4 * electrolyte solutions
5 * (see @ref thermoprops and @link Cantera::DebyeHuckel DebyeHuckel @endlink).
6 *
7 * Class DebyeHuckel represents a dilute liquid electrolyte phase which
8 * obeys the Debye Huckel formulation for nonideality.
9 */
10
11// This file is part of Cantera. See License.txt in the top-level directory or
12// at https://cantera.org/license.txt for license and copyright information.
13
21
22#include <cstdio>
23
24namespace Cantera
25{
26
27namespace {
28double A_Debye_default = 1.172576; // units = sqrt(kg/gmol)
29double B_Debye_default = 3.28640E9; // units = sqrt(kg/gmol) / m
30double maxIionicStrength_default = 30.0;
31}
32
33DebyeHuckel::DebyeHuckel(const string& inputFile, const string& id_)
34 : m_maxIionicStrength(maxIionicStrength_default)
35 , m_A_Debye(A_Debye_default)
36 , m_B_Debye(B_Debye_default)
37{
38 initThermoFile(inputFile, id_);
39}
40
41DebyeHuckel::~DebyeHuckel()
42{
43 // Defined in .cpp to limit dependence on WaterProps.h
44}
45
46// -------- Molar Thermodynamic Properties of the Solution ---------------
47
49{
51 return mean_X(m_tmpV);
52}
53
55{
57 return mean_X(m_tmpV);
58}
59
61{
63 return mean_X(m_tmpV);
64}
65
67{
69 return mean_X(m_tmpV);
70}
71
72// ------- Mechanical Equation of State Properties ------------------------
73
75{
76 if (m_waterSS) {
77 // Store the internal density of the water SS. Note, we would have to do
78 // this for all other species if they had pressure dependent properties.
80 }
82 double dd = meanMolecularWeight() / mean_X(m_tmpV);
84}
85
86// ------- Activities and Activity Concentrations
87
89{
90 double c_solvent = standardConcentration();
92 for (size_t k = 0; k < m_kk; k++) {
93 c[k] *= c_solvent;
94 }
95}
96
98{
99 double mvSolvent = providePDSS(0)->molarVolume();
100 return 1.0 / mvSolvent;
101}
102
103void DebyeHuckel::getActivities(double* ac) const
104{
106
107 // Update the molality array, m_molalities(). This requires an update due to
108 // mole fractions
110 for (size_t k = 1; k < m_kk; k++) {
111 ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
112 }
113 double xmolSolvent = moleFraction(0);
114 ac[0] = exp(m_lnActCoeffMolal[0]) * xmolSolvent;
115}
116
118{
120 A_Debye_TP(-1.0, -1.0);
122 copy(m_lnActCoeffMolal.begin(), m_lnActCoeffMolal.end(), acMolality);
123 for (size_t k = 0; k < m_kk; k++) {
124 acMolality[k] = exp(acMolality[k]);
125 }
126}
127
128// ------ Partial Molar Properties of the Solution -----------------
129
130void DebyeHuckel::getChemPotentials(double* mu) const
131{
132 double xx;
133
134 // First get the standard chemical potentials in molar form. This requires
135 // updates of standard state as a function of T and P
137
138 // Update the activity coefficients. This also updates the internal molality
139 // array.
141 double xmolSolvent = moleFraction(0);
142 for (size_t k = 1; k < m_kk; k++) {
143 xx = std::max(m_molalities[k], SmallNumber);
144 mu[k] += RT() * (log(xx) + m_lnActCoeffMolal[k]);
145 }
146 xx = std::max(xmolSolvent, SmallNumber);
147 mu[0] += RT() * (log(xx) + m_lnActCoeffMolal[0]);
148}
149
151{
152 // Get the nondimensional standard state enthalpies
153 getEnthalpy_RT(hbar);
154
155 // Dimensionalize it.
156 for (size_t k = 0; k < m_kk; k++) {
157 hbar[k] *= RT();
158 }
159
160 // Check to see whether activity coefficients are temperature
161 // dependent. If they are, then calculate the their temperature
162 // derivatives and add them into the result.
163 double dAdT = dA_DebyedT_TP();
164 if (dAdT != 0.0) {
165 // Update the activity coefficients, This also update the
166 // internally stored molalities.
169 for (size_t k = 0; k < m_kk; k++) {
170 hbar[k] -= RT() * temperature() * m_dlnActCoeffMolaldT[k];
171 }
172 }
173}
174
176{
177 // Get the standard state entropies at the temperature and pressure of the
178 // solution.
179 getEntropy_R(sbar);
180
181 // Dimensionalize the entropies
182 for (size_t k = 0; k < m_kk; k++) {
183 sbar[k] *= GasConstant;
184 }
185
186 // Update the activity coefficients, This also update the internally stored
187 // molalities.
189
190 // First we will add in the obvious dependence on the T term out front of
191 // the log activity term
192 double mm;
193 for (size_t k = 1; k < m_kk; k++) {
194 mm = std::max(SmallNumber, m_molalities[k]);
195 sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal[k]);
196 }
197 double xmolSolvent = moleFraction(0);
198 mm = std::max(SmallNumber, xmolSolvent);
199 sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal[0]);
200
201 // Check to see whether activity coefficients are temperature dependent. If
202 // they are, then calculate the their temperature derivatives and add them
203 // into the result.
204 double dAdT = dA_DebyedT_TP();
205 if (dAdT != 0.0) {
207 for (size_t k = 0; k < m_kk; k++) {
208 sbar[k] -= RT() * m_dlnActCoeffMolaldT[k];
209 }
210 }
211}
212
214{
215 getStandardVolumes(vbar);
216
217 // Update the derivatives wrt the activity coefficients.
220 for (size_t k = 0; k < m_kk; k++) {
221 vbar[k] += RT() * m_dlnActCoeffMolaldP[k];
222 }
223}
224
225void DebyeHuckel::getPartialMolarCp(double* cpbar) const
226{
227 getCp_R(cpbar);
228 for (size_t k = 0; k < m_kk; k++) {
229 cpbar[k] *= GasConstant;
230 }
231
232 // Check to see whether activity coefficients are temperature dependent. If
233 // they are, then calculate the their temperature derivatives and add them
234 // into the result.
235 double dAdT = dA_DebyedT_TP();
236 if (dAdT != 0.0) {
237 // Update the activity coefficients, This also update the internally
238 // stored molalities.
242 for (size_t k = 0; k < m_kk; k++) {
243 cpbar[k] -= (2.0 * RT() * m_dlnActCoeffMolaldT[k] +
245 }
246 }
247}
248
249// -------------- Utilities -------------------------------
250
251//! Utility function to assign an integer value from a string for the
252//! ElectrolyteSpeciesType field.
253/*!
254 * @param estString input string that will be interpreted
255 */
256static int interp_est(const string& estString)
257{
258 if (caseInsensitiveEquals(estString, "solvent")) {
259 return cEST_solvent;
260 } else if (estString == "charged-species"
261 || caseInsensitiveEquals(estString, "chargedspecies")) {
262 return cEST_chargedSpecies;
263 } else if (estString == "weak-acid-associated"
264 || caseInsensitiveEquals(estString, "weakacidassociated")) {
265 return cEST_weakAcidAssociated;
266 } else if (estString == "strong-acid-associated"
267 || caseInsensitiveEquals(estString, "strongacidassociated")) {
268 return cEST_strongAcidAssociated;
269 } else if (estString == "polar-neutral"
270 || caseInsensitiveEquals(estString, "polarneutral")) {
271 return cEST_polarNeutral;
272 } else if (estString == "nonpolar-neutral"
273 || caseInsensitiveEquals(estString, "nonpolarneutral")) {
274 return cEST_nonpolarNeutral;
275 } else {
276 throw CanteraError("interp_est (DebyeHuckel)",
277 "Invalid electrolyte species type '{}'", estString);
278 }
279}
280
281void DebyeHuckel::setDebyeHuckelModel(const string& model) {
282 if (model == ""
283 || model == "dilute-limit"
284 || caseInsensitiveEquals(model, "Dilute_limit")) {
285 m_formDH = DHFORM_DILUTE_LIMIT;
286 } else if (model == "B-dot-with-variable-a"
287 || caseInsensitiveEquals(model, "Bdot_with_variable_a")) {
288 m_formDH = DHFORM_BDOT_AK;
289 } else if (model == "B-dot-with-common-a"
290 || caseInsensitiveEquals(model, "Bdot_with_common_a")) {
291 m_formDH = DHFORM_BDOT_ACOMMON;
292 } else if (caseInsensitiveEquals(model, "beta_ij")) {
293 m_formDH = DHFORM_BETAIJ;
294 m_Beta_ij.resize(m_kk, m_kk, 0.0);
295 } else if (model == "Pitzer-with-beta_ij"
296 || caseInsensitiveEquals(model, "Pitzer_with_Beta_ij")) {
297 m_formDH = DHFORM_PITZER_BETAIJ;
298 m_Beta_ij.resize(m_kk, m_kk, 0.0);
299 } else {
300 throw CanteraError("DebyeHuckel::setDebyeHuckelModel",
301 "Unknown model '{}'", model);
302 }
303}
304
306{
307 if (A < 0) {
308 m_form_A_Debye = A_DEBYE_WATER;
309 } else {
310 m_form_A_Debye = A_DEBYE_CONST;
311 m_A_Debye = A;
312 }
313}
314
315void DebyeHuckel::setB_dot(double bdot)
316{
317 if (m_formDH == DHFORM_BETAIJ || m_formDH == DHFORM_DILUTE_LIMIT ||
318 m_formDH == DHFORM_PITZER_BETAIJ) {
319 throw CanteraError("DebyeHuckel::setB_dot",
320 "B_dot entry in the wrong DH form");
321 }
322 // Set B_dot parameters for charged species
323 for (size_t k = 0; k < nSpecies(); k++) {
324 if (fabs(charge(k)) > 0.0001) {
325 m_B_Dot[k] = bdot;
326 } else {
327 m_B_Dot[k] = 0.0;
328 }
329 }
330}
331
333{
334 m_Aionic_default = value;
335 for (size_t k = 0; k < m_kk; k++) {
336 if (std::isnan(m_Aionic[k])) {
337 m_Aionic[k] = value;
338 }
339 }
340}
341
342void DebyeHuckel::setBeta(const string& sp1, const string& sp2, double value)
343{
344 size_t k1 = speciesIndex(sp1);
345 if (k1 == npos) {
346 throw CanteraError("DebyeHuckel::setBeta", "Species '{}' not found", sp1);
347 }
348 size_t k2 = speciesIndex(sp2);
349 if (k2 == npos) {
350 throw CanteraError("DebyeHuckel::setBeta", "Species '{}' not found", sp2);
351 }
352 m_Beta_ij(k1, k2) = value;
353 m_Beta_ij(k2, k1) = value;
354}
355
357{
359 if (m_input.hasKey("activity-data")) {
360 auto& node = m_input["activity-data"].as<AnyMap>();
361 setDebyeHuckelModel(node["model"].asString());
362 if (node.hasKey("A_Debye")) {
363 if (node["A_Debye"] == "variable") {
364 setA_Debye(-1);
365 } else {
366 setA_Debye(node.convert("A_Debye", "kg^0.5/gmol^0.5"));
367 }
368 }
369 if (node.hasKey("B_Debye")) {
370 setB_Debye(node.convert("B_Debye", "kg^0.5/gmol^0.5/m"));
371 }
372 if (node.hasKey("max-ionic-strength")) {
373 setMaxIonicStrength(node["max-ionic-strength"].asDouble());
374 }
375 if (node.hasKey("use-Helgeson-fixed-form")) {
376 useHelgesonFixedForm(node["use-Helgeson-fixed-form"].asBool());
377 }
378 if (node.hasKey("default-ionic-radius")) {
379 setDefaultIonicRadius(node.convert("default-ionic-radius", "m"));
380 }
381 if (node.hasKey("B-dot")) {
382 setB_dot(node["B-dot"].asDouble());
383 }
384 if (node.hasKey("beta")) {
385 for (auto& item : node["beta"].asVector<AnyMap>()) {
386 auto& species = item["species"].asVector<string>(2);
387 setBeta(species[0], species[1], item["beta"].asDouble());
388 }
389 }
390 }
391
392 // Solvent
393 m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
394 if (m_waterSS) {
395 // Initialize the water property calculator. It will share the internal
396 // eos water calculator.
397 if (m_form_A_Debye == A_DEBYE_WATER) {
398 m_waterProps = make_unique<WaterProps>(m_waterSS);
399 }
400 } else if (dynamic_cast<PDSS_ConstVol*>(providePDSS(0)) == 0) {
401 throw CanteraError("DebyeHuckel::initThermo", "Solvent standard state"
402 " model must be WaterIAPWS or constant_incompressible.");
403 }
404
405 // Solutes
406 for (size_t k = 1; k < nSpecies(); k++) {
407 if (dynamic_cast<PDSS_ConstVol*>(providePDSS(k)) == 0) {
408 throw CanteraError("DebyeHuckel::initThermo", "Solute standard"
409 " state model must be constant_incompressible.");
410 }
411 }
412}
413
414void DebyeHuckel::getParameters(AnyMap& phaseNode) const
415{
417 AnyMap activityNode;
418
419 switch (m_formDH) {
420 case DHFORM_DILUTE_LIMIT:
421 activityNode["model"] = "dilute-limit";
422 break;
423 case DHFORM_BDOT_AK:
424 activityNode["model"] = "B-dot-with-variable-a";
425 break;
426 case DHFORM_BDOT_ACOMMON:
427 activityNode["model"] = "B-dot-with-common-a";
428 break;
429 case DHFORM_BETAIJ:
430 activityNode["model"] = "beta_ij";
431 break;
432 case DHFORM_PITZER_BETAIJ:
433 activityNode["model"] = "Pitzer-with-beta_ij";
434 break;
435 }
436
437 if (m_form_A_Debye == A_DEBYE_WATER) {
438 activityNode["A_Debye"] = "variable";
439 } else if (m_A_Debye != A_Debye_default) {
440 activityNode["A_Debye"].setQuantity(m_A_Debye, "kg^0.5/gmol^0.5");
441 }
442
443 if (m_B_Debye != B_Debye_default) {
444 activityNode["B_Debye"].setQuantity(m_B_Debye, "kg^0.5/gmol^0.5/m");
445 }
446 if (m_maxIionicStrength != maxIionicStrength_default) {
447 activityNode["max-ionic-strength"] = m_maxIionicStrength;
448 }
450 activityNode["use-Helgeson-fixed-form"] = true;
451 }
452 if (!isnan(m_Aionic_default)) {
453 activityNode["default-ionic-radius"].setQuantity(m_Aionic_default, "m");
454 }
455 for (double B_dot : m_B_Dot) {
456 if (B_dot != 0.0) {
457 activityNode["B-dot"] = B_dot;
458 break;
459 }
460 }
461 if (m_Beta_ij.nRows() && m_Beta_ij.nColumns()) {
462 vector<AnyMap> beta;
463 for (size_t i = 0; i < m_kk; i++) {
464 for (size_t j = i; j < m_kk; j++) {
465 if (m_Beta_ij(i, j) != 0) {
466 AnyMap entry;
467 entry["species"] = vector<string>{
468 speciesName(i), speciesName(j)};
469 entry["beta"] = m_Beta_ij(i, j);
470 beta.push_back(std::move(entry));
471 }
472 }
473 }
474 activityNode["beta"] = std::move(beta);
475 }
476 phaseNode["activity-data"] = std::move(activityNode);
477}
478
479void DebyeHuckel::getSpeciesParameters(const string& name, AnyMap& speciesNode) const
480{
482 size_t k = speciesIndex(name);
484 AnyMap dhNode;
485 if (m_Aionic[k] != m_Aionic_default) {
486 dhNode["ionic-radius"].setQuantity(m_Aionic[k], "m");
487 }
488
489 int estDefault = cEST_nonpolarNeutral;
490 if (k == 0) {
491 estDefault = cEST_solvent;
492 }
493
494 if (m_speciesCharge_Stoich[k] != charge(k)) {
495 dhNode["weak-acid-charge"] = m_speciesCharge_Stoich[k];
496 estDefault = cEST_weakAcidAssociated;
497 } else if (fabs(charge(k)) > 0.0001) {
498 estDefault = cEST_chargedSpecies;
499 }
500
501 if (m_electrolyteSpeciesType[k] != estDefault) {
502 string estType;
503 switch (m_electrolyteSpeciesType[k]) {
504 case cEST_solvent:
505 estType = "solvent";
506 break;
507 case cEST_chargedSpecies:
508 estType = "charged-species";
509 break;
510 case cEST_weakAcidAssociated:
511 estType = "weak-acid-associated";
512 break;
513 case cEST_strongAcidAssociated:
514 estType = "strong-acid-associated";
515 break;
516 case cEST_polarNeutral:
517 estType = "polar-neutral";
518 break;
519 case cEST_nonpolarNeutral:
520 estType = "nonpolar-neutral";
521 break;
522 default:
523 throw CanteraError("DebyeHuckel::getSpeciesParameters",
524 "Unknown electrolyte species type ({}) for species '{}'",
526 }
527 dhNode["electrolyte-species-type"] = estType;
528 }
529
530 if (dhNode.size()) {
531 speciesNode["Debye-Huckel"] = std::move(dhNode);
532 }
533}
534
535
536double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const
537{
538 double T = temperature();
539 double A;
540 if (tempArg != -1.0) {
541 T = tempArg;
542 }
543 double P = pressure();
544 if (presArg != -1.0) {
545 P = presArg;
546 }
547
548 switch (m_form_A_Debye) {
549 case A_DEBYE_CONST:
550 A = m_A_Debye;
551 break;
552 case A_DEBYE_WATER:
553 A = m_waterProps->ADebye(T, P, 0);
554 m_A_Debye = A;
555 break;
556 default:
557 throw CanteraError("DebyeHuckel::A_Debye_TP", "shouldn't be here");
558 }
559 return A;
560}
561
562double DebyeHuckel::dA_DebyedT_TP(double tempArg, double presArg) const
563{
564 double T = temperature();
565 if (tempArg != -1.0) {
566 T = tempArg;
567 }
568 double P = pressure();
569 if (presArg != -1.0) {
570 P = presArg;
571 }
572 double dAdT;
573 switch (m_form_A_Debye) {
574 case A_DEBYE_CONST:
575 dAdT = 0.0;
576 break;
577 case A_DEBYE_WATER:
578 dAdT = m_waterProps->ADebye(T, P, 1);
579 break;
580 default:
581 throw CanteraError("DebyeHuckel::dA_DebyedT_TP", "shouldn't be here");
582 }
583 return dAdT;
584}
585
586double DebyeHuckel::d2A_DebyedT2_TP(double tempArg, double presArg) const
587{
588 double T = temperature();
589 if (tempArg != -1.0) {
590 T = tempArg;
591 }
592 double P = pressure();
593 if (presArg != -1.0) {
594 P = presArg;
595 }
596 double d2AdT2;
597 switch (m_form_A_Debye) {
598 case A_DEBYE_CONST:
599 d2AdT2 = 0.0;
600 break;
601 case A_DEBYE_WATER:
602 d2AdT2 = m_waterProps->ADebye(T, P, 2);
603 break;
604 default:
605 throw CanteraError("DebyeHuckel::d2A_DebyedT2_TP", "shouldn't be here");
606 }
607 return d2AdT2;
608}
609
610double DebyeHuckel::dA_DebyedP_TP(double tempArg, double presArg) const
611{
612 double T = temperature();
613 if (tempArg != -1.0) {
614 T = tempArg;
615 }
616 double P = pressure();
617 if (presArg != -1.0) {
618 P = presArg;
619 }
620 double dAdP;
621 switch (m_form_A_Debye) {
622 case A_DEBYE_CONST:
623 dAdP = 0.0;
624 break;
625 case A_DEBYE_WATER:
626 dAdP = m_waterProps->ADebye(T, P, 3);
627 break;
628 default:
629 throw CanteraError("DebyeHuckel::dA_DebyedP_TP", "shouldn't be here");
630 }
631 return dAdP;
632}
633
634// ---------- Other Property Functions
635
636double DebyeHuckel::AionicRadius(int k) const
637{
638 return m_Aionic[k];
639}
640
641// ------------ Private and Restricted Functions ------------------
642
643bool DebyeHuckel::addSpecies(shared_ptr<Species> spec)
644{
645 bool added = MolalityVPSSTP::addSpecies(spec);
646 if (added) {
647 m_lnActCoeffMolal.push_back(0.0);
648 m_dlnActCoeffMolaldT.push_back(0.0);
649 m_d2lnActCoeffMolaldT2.push_back(0.0);
650 m_dlnActCoeffMolaldP.push_back(0.0);
651 m_B_Dot.push_back(0.0);
652 m_tmpV.push_back(0.0);
653
654 // NAN will be replaced with default value
655 double Aionic = NAN;
656
657 // Guess electrolyte species type based on charge properties
658 int est = cEST_nonpolarNeutral;
659 double stoichCharge = spec->charge;
660 if (fabs(spec->charge) > 0.0001) {
661 est = cEST_chargedSpecies;
662 }
663
664 if (spec->input.hasKey("Debye-Huckel")) {
665 auto& dhNode = spec->input["Debye-Huckel"].as<AnyMap>();
666 Aionic = dhNode.convert("ionic-radius", "m", NAN);
667 if (dhNode.hasKey("weak-acid-charge")) {
668 stoichCharge = dhNode["weak-acid-charge"].asDouble();
669 if (fabs(stoichCharge - spec->charge) > 0.0001) {
670 est = cEST_weakAcidAssociated;
671 }
672 }
673 // Apply override of the electrolyte species type
674 if (dhNode.hasKey("electrolyte-species-type")) {
675 est = interp_est(dhNode["electrolyte-species-type"].asString());
676 }
677 }
678
679 if (m_electrolyteSpeciesType.size() == 0) {
680 est = cEST_solvent; // species 0 is the solvent
681 }
682
683 m_Aionic.push_back(Aionic);
684 m_speciesCharge_Stoich.push_back(stoichCharge);
685 m_electrolyteSpeciesType.push_back(est);
686 }
687 return added;
688}
689
690double DebyeHuckel::_nonpolarActCoeff(double IionicMolality)
691{
692 // These are coefficients to describe the increase in activity coeff for
693 // non-polar molecules due to the electrolyte becoming stronger (the so-
694 // called salt-out effect)
695 const static double npActCoeff[] = {0.1127, -0.01049, 1.545E-3};
696 double I2 = IionicMolality * IionicMolality;
697 double l10actCoeff =
698 npActCoeff[0] * IionicMolality +
699 npActCoeff[1] * I2 +
700 npActCoeff[2] * I2 * IionicMolality;
701 return pow(10.0 , l10actCoeff);
702}
703
705{
706 const double a0 = 1.454;
707 const double b0 = 0.02236;
708 const double c0 = 9.380E-3;
709 const double d0 = -5.362E-4;
710 double Is = m_IionicMolalityStoich;
711 if (Is <= 0.0) {
712 return 0.0;
713 }
714 double Is2 = Is * Is;
715 double bhat = 1.0 + a0 * sqrt(Is);
716 double func = bhat - 2.0 * log(bhat) - 1.0/bhat;
717 double v1 = m_A_Debye / (a0 * a0 * a0 * Is) * func;
718 double oc = 1.0 - v1 + b0 * Is / 2.0 + 2.0 * c0 * Is2 / 3.0
719 + 3.0 * d0 * Is2 * Is / 4.0;
720 return oc;
721}
722
724{
725 // Update the internally stored vector of molalities
727 double oc = _osmoticCoeffHelgesonFixedForm();
728 double sum = 0.0;
729 for (size_t k = 1; k < m_kk; k++) {
730 sum += std::max(m_molalities[k], 0.0);
731 }
732 if (sum > 2.0 * m_maxIionicStrength) {
733 sum = 2.0 * m_maxIionicStrength;
734 };
735 return - m_Mnaught * sum * oc;
736}
737
739{
740 double z_k, zs_k1, zs_k2;
741
742 // Update the internally stored vector of molalities
744
745 // Calculate the apparent (real) ionic strength.
746 //
747 // Note this is not the stoichiometric ionic strengh, where reactions of
748 // ions forming neutral salts are ignorred in calculating the ionic
749 // strength.
750 m_IionicMolality = 0.0;
751 for (size_t k = 0; k < m_kk; k++) {
752 z_k = m_speciesCharge[k];
753 m_IionicMolality += m_molalities[k] * z_k * z_k;
754 }
755 m_IionicMolality /= 2.0;
757
758 // Calculate the stoichiometric ionic charge
760 for (size_t k = 0; k < m_kk; k++) {
761 z_k = m_speciesCharge[k];
762 zs_k1 = m_speciesCharge_Stoich[k];
763 if (z_k == zs_k1) {
764 m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
765 } else {
766 zs_k2 = z_k - zs_k1;
767 m_IionicMolalityStoich += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
768 }
769 }
772
773 // Possibly update the stored value of the Debye-Huckel parameter A_Debye
774 // This parameter appears on the top of the activity coefficient expression.
775 // It depends on T (and P), as it depends explicitly on the temperature.
776 // Also, the dielectric constant is usually a fairly strong function of T,
777 // also.
779
780 // Calculate a safe value for the mole fraction of the solvent
781 double xmolSolvent = moleFraction(0);
782 xmolSolvent = std::max(8.689E-3, xmolSolvent);
783
784 int est;
785 double ac_nonPolar = 1.0;
786 double numTmp = m_A_Debye * sqrt(m_IionicMolality);
787 double denomTmp = m_B_Debye * sqrt(m_IionicMolality);
788 double coeff;
789 double lnActivitySolvent = 0.0;
790 double tmp;
791 double tmpLn;
792 double y, yp1, sigma;
793 switch (m_formDH) {
794 case DHFORM_DILUTE_LIMIT:
795 for (size_t k = 0; k < m_kk; k++) {
796 z_k = m_speciesCharge[k];
797 m_lnActCoeffMolal[k] = - z_k * z_k * numTmp;
798 }
799 lnActivitySolvent =
800 (xmolSolvent - 1.0)/xmolSolvent +
801 2.0 / 3.0 * m_A_Debye * m_Mnaught *
803 break;
804
805 case DHFORM_BDOT_AK:
806 ac_nonPolar = _nonpolarActCoeff(m_IionicMolality);
807 for (size_t k = 0; k < m_kk; k++) {
809 if (est == cEST_nonpolarNeutral) {
810 m_lnActCoeffMolal[k] = log(ac_nonPolar);
811 } else {
812 z_k = m_speciesCharge[k];
814 - z_k * z_k * numTmp / (1.0 + denomTmp * m_Aionic[k])
815 + log(10.0) * m_B_Dot[k] * m_IionicMolality;
816 }
817 }
818
819 lnActivitySolvent = (xmolSolvent - 1.0)/xmolSolvent;
820 coeff = 2.0 / 3.0 * m_A_Debye * m_Mnaught
821 * sqrt(m_IionicMolality);
822 tmp = 0.0;
823 if (denomTmp > 0.0) {
824 for (size_t k = 0; k < m_kk; k++) {
825 if (k != 0 || m_Aionic[k] != 0.0) {
826 y = denomTmp * m_Aionic[k];
827 yp1 = y + 1.0;
828 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
829 z_k = m_speciesCharge[k];
830 tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
831 }
832 }
833 }
834 lnActivitySolvent += coeff * tmp;
835 tmp = 0.0;
836 for (size_t k = 1; k < m_kk; k++) {
837 z_k = m_speciesCharge[k];
838 if (z_k != 0.0) {
839 tmp += m_B_Dot[k] * m_molalities[k];
840 }
841 }
842 lnActivitySolvent -=
843 m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
844
845 // Special section to implement the Helgeson fixed form for the water
846 // brine activity coefficient.
848 lnActivitySolvent = _lnactivityWaterHelgesonFixedForm();
849 }
850 break;
851
852 case DHFORM_BDOT_ACOMMON:
853 denomTmp *= m_Aionic[0];
854 for (size_t k = 0; k < m_kk; k++) {
855 z_k = m_speciesCharge[k];
857 - z_k * z_k * numTmp / (1.0 + denomTmp)
858 + log(10.0) * m_B_Dot[k] * m_IionicMolality;
859 }
860 if (denomTmp > 0.0) {
861 y = denomTmp;
862 yp1 = y + 1.0;
863 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
864 } else {
865 sigma = 0.0;
866 }
867 lnActivitySolvent =
868 (xmolSolvent - 1.0)/xmolSolvent +
869 2.0 /3.0 * m_A_Debye * m_Mnaught *
870 m_IionicMolality * sqrt(m_IionicMolality) * sigma;
871 tmp = 0.0;
872 for (size_t k = 1; k < m_kk; k++) {
873 z_k = m_speciesCharge[k];
874 if (z_k != 0.0) {
875 tmp += m_B_Dot[k] * m_molalities[k];
876 }
877 }
878 lnActivitySolvent -=
879 m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
880 break;
881
882 case DHFORM_BETAIJ:
883 denomTmp = m_B_Debye * m_Aionic[0];
884 denomTmp *= sqrt(m_IionicMolality);
885 lnActivitySolvent =
886 (xmolSolvent - 1.0)/xmolSolvent;
887
888 for (size_t k = 1; k < m_kk; k++) {
889 z_k = m_speciesCharge[k];
891 - z_k * z_k * numTmp / (1.0 + denomTmp);
892 for (size_t j = 0; j < m_kk; j++) {
893 double beta = m_Beta_ij.value(k, j);
894 m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
895 }
896 }
897 if (denomTmp > 0.0) {
898 y = denomTmp;
899 yp1 = y + 1.0;
900 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1));
901 } else {
902 sigma = 0.0;
903 }
904 lnActivitySolvent =
905 (xmolSolvent - 1.0)/xmolSolvent +
906 2.0 /3.0 * m_A_Debye * m_Mnaught *
907 m_IionicMolality * sqrt(m_IionicMolality) * sigma;
908 tmp = 0.0;
909 for (size_t k = 0; k < m_kk; k++) {
910 for (size_t j = 0; j < m_kk; j++) {
911 tmp +=
913 }
914 }
915 lnActivitySolvent -= m_Mnaught * tmp;
916 break;
917
918 case DHFORM_PITZER_BETAIJ:
919 denomTmp = m_B_Debye * sqrt(m_IionicMolality);
920 denomTmp *= m_Aionic[0];
921 numTmp = m_A_Debye * sqrt(m_IionicMolality);
922 tmpLn = log(1.0 + denomTmp);
923 for (size_t k = 1; k < m_kk; k++) {
924 z_k = m_speciesCharge[k];
926 - z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
928 - 2.0 * z_k * z_k * m_A_Debye * tmpLn /
929 (3.0 * m_B_Debye * m_Aionic[0]);
930 for (size_t j = 0; j < m_kk; j++) {
931 m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
932 m_Beta_ij.value(k, j);
933 }
934 }
935 sigma = 1.0 / (1.0 + denomTmp);
936 lnActivitySolvent =
937 (xmolSolvent - 1.0)/xmolSolvent +
938 2.0 /3.0 * m_A_Debye * m_Mnaught *
939 m_IionicMolality * sqrt(m_IionicMolality) * sigma;
940 tmp = 0.0;
941 for (size_t k = 0; k < m_kk; k++) {
942 for (size_t j = 0; j < m_kk; j++) {
943 tmp +=
945 }
946 }
947 lnActivitySolvent -= m_Mnaught * tmp;
948 break;
949
950 default:
951 throw CanteraError("DebyeHuckel::s_update_lnMolalityActCoeff", "ERROR");
952 }
953
954 // Above, we calculated the ln(activitySolvent). Translate that into the
955 // molar-based activity coefficient by dividing by the solvent mole
956 // fraction. Solvents are not on the molality scale.
957 xmolSolvent = moleFraction(0);
958 m_lnActCoeffMolal[0] = lnActivitySolvent - log(xmolSolvent);
959}
960
962{
963 double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
964 // First we store dAdT explicitly here
965 double dAdT = dA_DebyedT_TP();
966 if (dAdT == 0.0) {
967 for (size_t k = 0; k < m_kk; k++) {
968 m_dlnActCoeffMolaldT[k] = 0.0;
969 }
970 return;
971 }
972
973 // Calculate a safe value for the mole fraction of the solvent
974 double xmolSolvent = moleFraction(0);
975 xmolSolvent = std::max(8.689E-3, xmolSolvent);
976 double sqrtI = sqrt(m_IionicMolality);
977 double numdAdTTmp = dAdT * sqrtI;
978 double denomTmp = m_B_Debye * sqrtI;
979 double d_lnActivitySolvent_dT = 0;
980
981 switch (m_formDH) {
982 case DHFORM_DILUTE_LIMIT:
983 for (size_t k = 1; k < m_kk; k++) {
985 m_lnActCoeffMolal[k] * dAdT / m_A_Debye;
986 }
987 d_lnActivitySolvent_dT = 2.0 / 3.0 * dAdT * m_Mnaught *
989 m_dlnActCoeffMolaldT[0] = d_lnActivitySolvent_dT;
990 break;
991
992 case DHFORM_BDOT_AK:
993 for (size_t k = 0; k < m_kk; k++) {
994 z_k = m_speciesCharge[k];
996 - z_k * z_k * numdAdTTmp / (1.0 + denomTmp * m_Aionic[k]);
997 }
998
999 m_dlnActCoeffMolaldT[0] = 0.0;
1000 coeff = 2.0 / 3.0 * dAdT * m_Mnaught * sqrtI;
1001 tmp = 0.0;
1002 if (denomTmp > 0.0) {
1003 for (size_t k = 0; k < m_kk; k++) {
1004 y = denomTmp * m_Aionic[k];
1005 yp1 = y + 1.0;
1006 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1007 z_k = m_speciesCharge[k];
1008 tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
1009 }
1010 }
1011 m_dlnActCoeffMolaldT[0] += coeff * tmp;
1012 break;
1013
1014 case DHFORM_BDOT_ACOMMON:
1015 denomTmp *= m_Aionic[0];
1016 for (size_t k = 0; k < m_kk; k++) {
1017 z_k = m_speciesCharge[k];
1019 - z_k * z_k * numdAdTTmp / (1.0 + denomTmp);
1020 }
1021 if (denomTmp > 0.0) {
1022 y = denomTmp;
1023 yp1 = y + 1.0;
1024 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1025 } else {
1026 sigma = 0.0;
1027 }
1028 m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
1029 m_IionicMolality * sqrtI * sigma;
1030 break;
1031
1032 case DHFORM_BETAIJ:
1033 denomTmp *= m_Aionic[0];
1034 for (size_t k = 1; k < m_kk; k++) {
1035 z_k = m_speciesCharge[k];
1036 m_dlnActCoeffMolaldT[k] = -z_k*z_k * numdAdTTmp / (1.0 + denomTmp);
1037 }
1038 if (denomTmp > 0.0) {
1039 y = denomTmp;
1040 yp1 = y + 1.0;
1041 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1042 } else {
1043 sigma = 0.0;
1044 }
1045 m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
1046 m_IionicMolality * sqrtI * sigma;
1047 break;
1048
1049 case DHFORM_PITZER_BETAIJ:
1050 denomTmp *= m_Aionic[0];
1051 tmpLn = log(1.0 + denomTmp);
1052 for (size_t k = 1; k < m_kk; k++) {
1053 z_k = m_speciesCharge[k];
1055 - z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
1056 - 2.0 * z_k * z_k * dAdT * tmpLn / (m_B_Debye * m_Aionic[0]);
1057 m_dlnActCoeffMolaldT[k] /= 3.0;
1058 }
1059
1060 sigma = 1.0 / (1.0 + denomTmp);
1061 m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
1062 m_IionicMolality * sqrtI * sigma;
1063 break;
1064
1065 default:
1066 throw CanteraError("DebyeHuckel::s_update_dlnMolalityActCoeff_dT",
1067 "ERROR");
1068 }
1069}
1070
1072{
1073 double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
1074 double dAdT = dA_DebyedT_TP();
1075 double d2AdT2 = d2A_DebyedT2_TP();
1076 if (d2AdT2 == 0.0 && dAdT == 0.0) {
1077 for (size_t k = 0; k < m_kk; k++) {
1078 m_d2lnActCoeffMolaldT2[k] = 0.0;
1079 }
1080 return;
1081 }
1082
1083 // Calculate a safe value for the mole fraction of the solvent
1084 double xmolSolvent = moleFraction(0);
1085 xmolSolvent = std::max(8.689E-3, xmolSolvent);
1086 double sqrtI = sqrt(m_IionicMolality);
1087 double numd2AdT2Tmp = d2AdT2 * sqrtI;
1088 double denomTmp = m_B_Debye * sqrtI;
1089
1090 switch (m_formDH) {
1091 case DHFORM_DILUTE_LIMIT:
1092 for (size_t k = 0; k < m_kk; k++) {
1094 m_lnActCoeffMolal[k] * d2AdT2 / m_A_Debye;
1095 }
1096 break;
1097
1098 case DHFORM_BDOT_AK:
1099 for (size_t k = 0; k < m_kk; k++) {
1100 z_k = m_speciesCharge[k];
1102 - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp * m_Aionic[k]);
1103 }
1104
1105 m_d2lnActCoeffMolaldT2[0] = 0.0;
1106 coeff = 2.0 / 3.0 * d2AdT2 * m_Mnaught * sqrtI;
1107 tmp = 0.0;
1108 if (denomTmp > 0.0) {
1109 for (size_t k = 0; k < m_kk; k++) {
1110 y = denomTmp * m_Aionic[k];
1111 yp1 = y + 1.0;
1112 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1113 z_k = m_speciesCharge[k];
1114 tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
1115 }
1116 }
1117 m_d2lnActCoeffMolaldT2[0] += coeff * tmp;
1118 break;
1119
1120 case DHFORM_BDOT_ACOMMON:
1121 denomTmp *= m_Aionic[0];
1122 for (size_t k = 0; k < m_kk; k++) {
1123 z_k = m_speciesCharge[k];
1125 - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp);
1126 }
1127 if (denomTmp > 0.0) {
1128 y = denomTmp;
1129 yp1 = y + 1.0;
1130 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1131 } else {
1132 sigma = 0.0;
1133 }
1134 m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
1135 m_IionicMolality * sqrtI * sigma;
1136 break;
1137
1138 case DHFORM_BETAIJ:
1139 denomTmp *= m_Aionic[0];
1140 for (size_t k = 1; k < m_kk; k++) {
1141 z_k = m_speciesCharge[k];
1142 m_d2lnActCoeffMolaldT2[k] = -z_k*z_k * numd2AdT2Tmp / (1.0 + denomTmp);
1143 }
1144 if (denomTmp > 0.0) {
1145 y = denomTmp;
1146 yp1 = y + 1.0;
1147 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1));
1148 } else {
1149 sigma = 0.0;
1150 }
1151 m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
1152 m_IionicMolality * sqrtI * sigma;
1153 break;
1154
1155 case DHFORM_PITZER_BETAIJ:
1156 denomTmp *= m_Aionic[0];
1157 tmpLn = log(1.0 + denomTmp);
1158 for (size_t k = 1; k < m_kk; k++) {
1159 z_k = m_speciesCharge[k];
1161 - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
1162 - 2.0 * z_k * z_k * d2AdT2 * tmpLn / (m_B_Debye * m_Aionic[0]);
1163 m_d2lnActCoeffMolaldT2[k] /= 3.0;
1164 }
1165
1166 sigma = 1.0 / (1.0 + denomTmp);
1167 m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
1168 m_IionicMolality * sqrtI * sigma;
1169 break;
1170
1171 default:
1172 throw CanteraError("DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2",
1173 "ERROR");
1174 }
1175}
1176
1178{
1179 double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
1180 int est;
1181 double dAdP = dA_DebyedP_TP();
1182 if (dAdP == 0.0) {
1183 for (size_t k = 0; k < m_kk; k++) {
1184 m_dlnActCoeffMolaldP[k] = 0.0;
1185 }
1186 return;
1187 }
1188
1189 // Calculate a safe value for the mole fraction of the solvent
1190 double xmolSolvent = moleFraction(0);
1191 xmolSolvent = std::max(8.689E-3, xmolSolvent);
1192 double sqrtI = sqrt(m_IionicMolality);
1193 double numdAdPTmp = dAdP * sqrtI;
1194 double denomTmp = m_B_Debye * sqrtI;
1195
1196 switch (m_formDH) {
1197 case DHFORM_DILUTE_LIMIT:
1198 for (size_t k = 0; k < m_kk; k++) {
1200 m_lnActCoeffMolal[k] * dAdP / m_A_Debye;
1201 }
1202 break;
1203
1204 case DHFORM_BDOT_AK:
1205 for (size_t k = 0; k < m_kk; k++) {
1206 est = m_electrolyteSpeciesType[k];
1207 if (est == cEST_nonpolarNeutral) {
1208 m_lnActCoeffMolal[k] = 0.0;
1209 } else {
1210 z_k = m_speciesCharge[k];
1212 - z_k * z_k * numdAdPTmp / (1.0 + denomTmp * m_Aionic[k]);
1213 }
1214 }
1215
1216 m_dlnActCoeffMolaldP[0] = 0.0;
1217 coeff = 2.0 / 3.0 * dAdP * m_Mnaught * sqrtI;
1218 tmp = 0.0;
1219 if (denomTmp > 0.0) {
1220 for (size_t k = 0; k < m_kk; k++) {
1221 y = denomTmp * m_Aionic[k];
1222 yp1 = y + 1.0;
1223 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1224 z_k = m_speciesCharge[k];
1225 tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
1226 }
1227 }
1228 m_dlnActCoeffMolaldP[0] += coeff * tmp;
1229 break;
1230
1231 case DHFORM_BDOT_ACOMMON:
1232 denomTmp *= m_Aionic[0];
1233 for (size_t k = 0; k < m_kk; k++) {
1234 z_k = m_speciesCharge[k];
1236 - z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
1237 }
1238 if (denomTmp > 0.0) {
1239 y = denomTmp;
1240 yp1 = y + 1.0;
1241 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1242 } else {
1243 sigma = 0.0;
1244 }
1246 2.0 /3.0 * dAdP * m_Mnaught *
1247 m_IionicMolality * sqrtI * sigma;
1248 break;
1249
1250 case DHFORM_BETAIJ:
1251 denomTmp *= m_Aionic[0];
1252 for (size_t k = 1; k < m_kk; k++) {
1253 z_k = m_speciesCharge[k];
1254 m_dlnActCoeffMolaldP[k] = - z_k*z_k * numdAdPTmp / (1.0 + denomTmp);
1255 }
1256 if (denomTmp > 0.0) {
1257 y = denomTmp;
1258 yp1 = y + 1.0;
1259 sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1260 } else {
1261 sigma = 0.0;
1262 }
1263 m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
1264 m_IionicMolality * sqrtI * sigma;
1265 break;
1266
1267 case DHFORM_PITZER_BETAIJ:
1268 denomTmp *= m_Aionic[0];
1269 tmpLn = log(1.0 + denomTmp);
1270 for (size_t k = 1; k < m_kk; k++) {
1271 z_k = m_speciesCharge[k];
1273 - z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
1274 - 2.0 * z_k * z_k * dAdP * tmpLn
1275 / (m_B_Debye * m_Aionic[0]);
1276 m_dlnActCoeffMolaldP[k] /= 3.0;
1277 }
1278
1279 sigma = 1.0 / (1.0 + denomTmp);
1280 m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
1281 m_IionicMolality * sqrtI * sigma;
1282 break;
1283
1284 default:
1285 throw CanteraError("DebyeHuckel::s_update_dlnMolalityActCoeff_dP",
1286 "ERROR");
1287 }
1288}
1289
1290}
Headers for the DebyeHuckel ThermoPhase object, which models dilute electrolyte solutions (see Thermo...
Declarations for the class PDSS_ConstVol (pressure dependent standard state) which handles calculatio...
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
Declaration for class Cantera::Species.
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 convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
Definition AnyMap.cpp:1535
size_t nRows() const
Number of rows.
Definition Array.h:176
size_t nColumns() const
Number of columns.
Definition Array.h:181
double & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Definition Array.h:160
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.
Array2D m_Beta_ij
Array of 2D data used in the DHFORM_BETAIJ formulation Beta_ij.value(i,j) is the coefficient of the j...
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
int m_formDH
form of the Debye-Huckel parameterization used in the model.
DebyeHuckel(const string &inputFile="", const string &id="")
Full constructor for creating the phase.
double m_A_Debye
Current value of the Debye Constant, A_Debye.
virtual double d2A_DebyedT2_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the 2nd derivative of the Debye Huckel constant with respect to temperature as a function of...
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
vector< double > m_B_Dot
Array of B_Dot values.
double m_IionicMolality
Current value of the ionic strength on the molality scale.
double _osmoticCoeffHelgesonFixedForm() const
Formula for the osmotic coefficient that occurs in the GWB.
double m_densWaterSS
Storage for the density of water's standard state.
void s_update_d2lnMolalityActCoeff_dT2() const
Calculate the temperature 2nd derivative of the activity coefficient.
void getSpeciesParameters(const string &name, AnyMap &speciesNode) const override
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
bool m_useHelgesonFixedForm
If true, then the fixed for of Helgeson's activity for water is used instead of the rigorous form obt...
static double _nonpolarActCoeff(double IionicMolality)
Static function that implements the non-polar species salt-out modifications.
vector< int > m_electrolyteSpeciesType
Vector containing the electrolyte species type.
double m_B_Debye
Current value of the constant that appears in the denominator.
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
double AionicRadius(int k=0) const
Reports the ionic radius of the kth species.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
vector< double > m_tmpV
vector of size m_kk, used as a temporary holding area.
void getActivityConcentrations(double *c) const override
This method returns an array of generalized concentrations.
void s_update_dlnMolalityActCoeff_dT() const
Calculation of temperature derivative of activity coefficient.
void s_update_lnMolalityActCoeff() const
Calculate the log activity coefficients.
vector< double > m_Aionic
a_k = Size of the ionic species in the DH formulation. units = meters
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
void setA_Debye(double A)
Set the A_Debye parameter.
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
void calcDensity() override
Calculate the density of the mixture using the partial molar volumes and mole fractions as input.
vector< double > m_dlnActCoeffMolaldT
Derivative of log act coeff wrt T.
double m_Aionic_default
Default ionic radius for species where it is not specified.
void setDebyeHuckelModel(const string &form)
Set the DebyeHuckel parameterization form.
vector< double > m_speciesCharge_Stoich
Stoichiometric species charge -> This is for calculations of the ionic strength which ignore ion-ion ...
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
double _lnactivityWaterHelgesonFixedForm() const
Formula for the log of the water activity that occurs in the GWB.
void setBeta(const string &sp1, const string &sp2, double value)
Set the value for the beta interaction between species sp1 and sp2.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
void getActivities(double *ac) const override
Get the array of non-dimensional activities at the current solution temperature, pressure,...
virtual double dA_DebyedT_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to temperature.
void getPartialMolarCp(double *cpbar) const override
Return an array of partial molar heat capacities for the species in the mixture.
void setDefaultIonicRadius(double value)
Set the default ionic radius [m] for each species.
vector< double > m_d2lnActCoeffMolaldT2
2nd Derivative of log act coeff wrt T
double gibbs_mole() const override
Molar Gibbs function. Units: J/kmol.
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
void s_update_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the activity coefficient.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
PDSS_Water * m_waterSS
Pointer to the Water standard state object.
double m_IionicMolalityStoich
Stoichiometric ionic strength on the molality scale.
vector< double > m_lnActCoeffMolal
Logarithm of the activity coefficients on the molality scale.
vector< double > m_dlnActCoeffMolaldP
Derivative of log act coeff wrt P.
void getMolalityActivityCoefficients(double *acMolality) const override
Get the array of non-dimensional molality-based activity coefficients at the current solution tempera...
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
virtual double A_Debye_TP(double temperature=-1.0, double pressure=-1.0) const
Return the Debye Huckel constant as a function of temperature and pressure (Units = sqrt(kg/gmol))
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...
double m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
vector< double > m_molalities
Current value of the molalities of the species in the phase.
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
Class for pressure dependent standard states that use a constant volume model.
Class for the liquid water pressure dependent standard state.
Definition PDSS_Water.h:50
double density() const override
Return the standard state density at standard state.
virtual double molarVolume() const
Return the molar volume at standard state.
Definition PDSS.cpp:63
void assignDensity(const double density_)
Set the internally stored constant density (kg/m^3) of the phase.
Definition Phase.cpp:597
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition Phase.cpp:153
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
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
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
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
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
vector< double > m_speciesCharge
Vector of species charges. length m_kk.
Definition Phase.h:853
string name() const
Return the name of the phase.
Definition Phase.cpp:20
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.
void initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
virtual void getSpeciesParameters(const string &name, AnyMap &speciesNode) const
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
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...
virtual void _updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
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.
Header file for a common definitions used in electrolytes thermodynamics.
bool caseInsensitiveEquals(const string &input, const string &test)
Case insensitive equality predicate.
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
static int interp_est(const string &estString)
Utility function to assign an integer value from a string for the ElectrolyteSpeciesType field.
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
const int cEST_solvent
Electrolyte species type.
Contains declarations for string manipulation functions within Cantera.