Cantera  3.1.0a1
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 
16 #include "cantera/thermo/Species.h"
21 
22 #include <cstdio>
23 
24 namespace Cantera
25 {
26 
27 namespace {
28 double A_Debye_default = 1.172576; // units = sqrt(kg/gmol)
29 double B_Debye_default = 3.28640E9; // units = sqrt(kg/gmol) / m
30 double maxIionicStrength_default = 30.0;
31 }
32 
33 DebyeHuckel::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 
41 DebyeHuckel::~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 {
62  getChemPotentials(m_tmpV.data());
63  return mean_X(m_tmpV);
64 }
65 
66 double DebyeHuckel::cp_mole() const
67 {
68  getPartialMolarCp(m_tmpV.data());
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();
91  getActivities(c);
92  for (size_t k = 0; k < m_kk; k++) {
93  c[k] *= c_solvent;
94  }
95 }
96 
97 double DebyeHuckel::standardConcentration(size_t k) const
98 {
99  double mvSolvent = providePDSS(0)->molarVolume();
100  return 1.0 / mvSolvent;
101 }
102 
103 void 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 
117 void DebyeHuckel::getMolalityActivityCoefficients(double* acMolality) const
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 
130 void 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 
213 void DebyeHuckel::getPartialMolarVolumes(double* vbar) const
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 
225 void 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  */
256 static 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 
281 void 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 
315 void 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 
342 void 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 
414 void 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 
479 void 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 
536 double 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 
562 double 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 
586 double 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 
610 double 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 
636 double DebyeHuckel::AionicRadius(int k) const
637 {
638  return m_Aionic[k];
639 }
640 
641 // ------------ Private and Restricted Functions ------------------
642 
643 bool 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 
690 double 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
726  calcMolalities();
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
743  calcMolalities();
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  }
770  m_IionicMolalityStoich /= 2.0;
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.
778  m_A_Debye = A_Debye_TP();
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++) {
808  est = m_electrolyteSpeciesType[k];
809  if (est == cEST_nonpolarNeutral) {
810  m_lnActCoeffMolal[k] = log(ac_nonPolar);
811  } else {
812  z_k = m_speciesCharge[k];
813  m_lnActCoeffMolal[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];
856  m_lnActCoeffMolal[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];
890  m_lnActCoeffMolal[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 +=
912  m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
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];
925  m_lnActCoeffMolal[k] =
926  - z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
927  m_lnActCoeffMolal[k] +=
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 +=
944  m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
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
double & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Definition: Array.h:160
size_t nRows() const
Number of rows.
Definition: Array.h:176
size_t nColumns() const
Number of columns.
Definition: Array.h:181
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.
Definition: ctexceptions.h:66
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...
Definition: DebyeHuckel.h:937
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
Definition: DebyeHuckel.cpp:48
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition: DebyeHuckel.h:912
int m_formDH
form of the Debye-Huckel parameterization used in the model.
Definition: DebyeHuckel.h:794
DebyeHuckel(const string &inputFile="", const string &id="")
Full constructor for creating the phase.
Definition: DebyeHuckel.cpp:33
double m_A_Debye
Current value of the Debye Constant, A_Debye.
Definition: DebyeHuckel.h:871
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.
Definition: DebyeHuckel.h:897
double m_IionicMolality
Current value of the ionic strength on the molality scale.
Definition: DebyeHuckel.h:816
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.
Definition: DebyeHuckel.h:909
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.
Definition: DebyeHuckel.h:849
bool m_useHelgesonFixedForm
If true, then the fixed for of Helgeson's activity for water is used instead of the rigorous form obt...
Definition: DebyeHuckel.h:827
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.
Definition: DebyeHuckel.h:807
double m_B_Debye
Current value of the constant that appears in the denominator.
Definition: DebyeHuckel.h:889
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.
Definition: DebyeHuckel.h:915
void getActivityConcentrations(double *c) const override
This method returns an array of generalized concentrations.
Definition: DebyeHuckel.cpp:88
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
Definition: DebyeHuckel.h:810
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.
Definition: DebyeHuckel.cpp:54
void calcDensity() override
Calculate the density of the mixture using the partial molar volumes and mole fractions as input.
Definition: DebyeHuckel.cpp:74
vector< double > m_dlnActCoeffMolaldT
Derivative of log act coeff wrt T.
Definition: DebyeHuckel.h:947
double m_Aionic_default
Default ionic radius for species where it is not specified.
Definition: DebyeHuckel.h:813
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 ...
Definition: DebyeHuckel.h:929
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition: DebyeHuckel.h:820
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.
Definition: DebyeHuckel.cpp:66
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
Definition: DebyeHuckel.h:950
double gibbs_mole() const override
Molar Gibbs function. Units: J/kmol.
Definition: DebyeHuckel.cpp:60
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
Definition: DebyeHuckel.cpp:97
void s_update_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the activity coefficient.
PDSS_Water * m_waterSS
Pointer to the Water standard state object.
Definition: DebyeHuckel.h:903
double m_IionicMolalityStoich
Stoichiometric ionic strength on the molality scale.
Definition: DebyeHuckel.h:830
vector< double > m_lnActCoeffMolal
Logarithm of the activity coefficients on the molality scale.
Definition: DebyeHuckel.h:944
vector< double > m_dlnActCoeffMolaldP
Derivative of log act coeff wrt P.
Definition: DebyeHuckel.h:953
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.
Class for pressure dependent standard states that use a constant volume model.
Definition: PDSS_ConstVol.h:23
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.
Definition: PDSS_Water.cpp:207
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.
Definition: ThermoPhase.h:1062
void initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
AnyMap m_input
Data supplied via setParameters.
Definition: ThermoPhase.h:1966
virtual bool addSpecies(shared_ptr< Species > spec)
Add a Species to this Phase.
Definition: Phase.cpp:701
double pressure() const override
Returns the current pressure of the phase.
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...
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
static int interp_est(const string &estString)
Utility function to assign an integer value from a string for the ElectrolyteSpeciesType field.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:158
const int cEST_solvent
Electrolyte species type.
Definition: electrolytes.h:18
Contains declarations for string manipulation functions within Cantera.