Cantera  2.4.0
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 http://www.cantera.org/license.txt for license and copyright information.
13 
20 #include "cantera/base/ctml.h"
21 
22 #include <cstdio>
23 
24 using namespace std;
25 
26 namespace Cantera
27 {
28 
29 DebyeHuckel::DebyeHuckel() :
30  m_formDH(DHFORM_DILUTE_LIMIT),
31  m_IionicMolality(0.0),
32  m_maxIionicStrength(30.0),
33  m_useHelgesonFixedForm(false),
34  m_IionicMolalityStoich(0.0),
35  m_form_A_Debye(A_DEBYE_CONST),
36  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
37  m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
38  m_waterSS(0),
39  m_densWaterSS(1000.)
40 {
41 }
42 
43 DebyeHuckel::DebyeHuckel(const std::string& inputFile,
44  const std::string& id_) :
45  m_formDH(DHFORM_DILUTE_LIMIT),
46  m_IionicMolality(0.0),
47  m_maxIionicStrength(30.0),
48  m_useHelgesonFixedForm(false),
49  m_IionicMolalityStoich(0.0),
50  m_form_A_Debye(A_DEBYE_CONST),
51  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
52  m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
53  m_waterSS(0),
54  m_densWaterSS(1000.)
55 {
56  initThermoFile(inputFile, id_);
57 }
58 
59 DebyeHuckel::DebyeHuckel(XML_Node& phaseRoot, const std::string& id_) :
60  m_formDH(DHFORM_DILUTE_LIMIT),
61  m_IionicMolality(0.0),
62  m_maxIionicStrength(3.0),
63  m_useHelgesonFixedForm(false),
64  m_IionicMolalityStoich(0.0),
65  m_form_A_Debye(A_DEBYE_CONST),
66  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
67  m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
68  m_waterSS(0),
69  m_densWaterSS(1000.)
70 {
71  importPhase(phaseRoot, this);
72 }
73 
74 DebyeHuckel::~DebyeHuckel()
75 {
76 }
77 
78 // -------- Molar Thermodynamic Properties of the Solution ---------------
79 
80 doublereal DebyeHuckel::enthalpy_mole() const
81 {
83  return mean_X(m_tmpV);
84 }
85 
86 doublereal DebyeHuckel::entropy_mole() const
87 {
89  return mean_X(m_tmpV);
90 }
91 
92 doublereal DebyeHuckel::gibbs_mole() const
93 {
94  getChemPotentials(m_tmpV.data());
95  return mean_X(m_tmpV);
96 }
97 
98 doublereal DebyeHuckel::cp_mole() const
99 {
100  getPartialMolarCp(m_tmpV.data());
101  return mean_X(m_tmpV);
102 }
103 
104 // ------- Mechanical Equation of State Properties ------------------------
105 
107 {
108  if (m_waterSS) {
109  // Store the internal density of the water SS. Note, we would have to do
110  // this for all other species if they had pressure dependent properties.
112  }
114  double dd = meanMolecularWeight() / mean_X(m_tmpV);
115  Phase::setDensity(dd);
116 }
117 
118 void DebyeHuckel::setDensity(doublereal rho)
119 {
120  double dens = density();
121  if (rho != dens) {
122  throw CanteraError("Idea;MolalSoln::setDensity",
123  "Density is not an independent variable");
124  }
125 }
126 
127 void DebyeHuckel::setMolarDensity(const doublereal conc)
128 {
129  double concI = molarDensity();
130  if (conc != concI) {
131  throw CanteraError("Idea;MolalSoln::setMolarDensity",
132  "molarDensity/density is not an independent variable");
133  }
134 }
135 
136 // ------- Activities and Activity Concentrations
137 
138 void DebyeHuckel::getActivityConcentrations(doublereal* c) const
139 {
140  double c_solvent = standardConcentration();
141  getActivities(c);
142  for (size_t k = 0; k < m_kk; k++) {
143  c[k] *= c_solvent;
144  }
145 }
146 
147 doublereal DebyeHuckel::standardConcentration(size_t k) const
148 {
149  double mvSolvent = providePDSS(0)->molarVolume();
150  return 1.0 / mvSolvent;
151 }
152 
153 void DebyeHuckel::getActivities(doublereal* ac) const
154 {
156 
157  // Update the molality array, m_molalities(). This requires an update due to
158  // mole fractions
160  for (size_t k = 1; k < m_kk; k++) {
161  ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
162  }
163  double xmolSolvent = moleFraction(0);
164  ac[0] = exp(m_lnActCoeffMolal[0]) * xmolSolvent;
165 }
166 
167 void DebyeHuckel::getMolalityActivityCoefficients(doublereal* acMolality) const
168 {
170  A_Debye_TP(-1.0, -1.0);
172  copy(m_lnActCoeffMolal.begin(), m_lnActCoeffMolal.end(), acMolality);
173  for (size_t k = 0; k < m_kk; k++) {
174  acMolality[k] = exp(acMolality[k]);
175  }
176 }
177 
178 // ------ Partial Molar Properties of the Solution -----------------
179 
180 void DebyeHuckel::getChemPotentials(doublereal* mu) const
181 {
182  double xx;
183 
184  // First get the standard chemical potentials in molar form. This requires
185  // updates of standard state as a function of T and P
187 
188  // Update the activity coefficients. This also updates the internal molality
189  // array.
191  double xmolSolvent = moleFraction(0);
192  for (size_t k = 1; k < m_kk; k++) {
193  xx = std::max(m_molalities[k], SmallNumber);
194  mu[k] += RT() * (log(xx) + m_lnActCoeffMolal[k]);
195  }
196  xx = std::max(xmolSolvent, SmallNumber);
197  mu[0] += RT() * (log(xx) + m_lnActCoeffMolal[0]);
198 }
199 
200 void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const
201 {
202  // Get the nondimensional standard state enthalpies
203  getEnthalpy_RT(hbar);
204 
205  // Dimensionalize it.
206  for (size_t k = 0; k < m_kk; k++) {
207  hbar[k] *= RT();
208  }
209 
210  // Check to see whether activity coefficients are temperature
211  // dependent. If they are, then calculate the their temperature
212  // derivatives and add them into the result.
213  double dAdT = dA_DebyedT_TP();
214  if (dAdT != 0.0) {
215  // Update the activity coefficients, This also update the
216  // internally stored molalities.
219  for (size_t k = 0; k < m_kk; k++) {
220  hbar[k] -= RT() * temperature() * m_dlnActCoeffMolaldT[k];
221  }
222  }
223 }
224 
225 void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const
226 {
227  // Get the standard state entropies at the temperature and pressure of the
228  // solution.
229  getEntropy_R(sbar);
230 
231  // Dimensionalize the entropies
232  for (size_t k = 0; k < m_kk; k++) {
233  sbar[k] *= GasConstant;
234  }
235 
236  // Update the activity coefficients, This also update the internally stored
237  // molalities.
239 
240  // First we will add in the obvious dependence on the T term out front of
241  // the log activity term
242  doublereal mm;
243  for (size_t k = 1; k < m_kk; k++) {
244  mm = std::max(SmallNumber, m_molalities[k]);
245  sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal[k]);
246  }
247  double xmolSolvent = moleFraction(0);
248  mm = std::max(SmallNumber, xmolSolvent);
249  sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal[0]);
250 
251  // Check to see whether activity coefficients are temperature dependent. If
252  // they are, then calculate the their temperature derivatives and add them
253  // into the result.
254  double dAdT = dA_DebyedT_TP();
255  if (dAdT != 0.0) {
257  for (size_t k = 0; k < m_kk; k++) {
258  sbar[k] -= RT() * m_dlnActCoeffMolaldT[k];
259  }
260  }
261 }
262 
263 void DebyeHuckel::getPartialMolarVolumes(doublereal* vbar) const
264 {
265  getStandardVolumes(vbar);
266 
267  // Update the derivatives wrt the activity coefficients.
270  for (size_t k = 0; k < m_kk; k++) {
271  vbar[k] += RT() * m_dlnActCoeffMolaldP[k];
272  }
273 }
274 
275 void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const
276 {
277  getCp_R(cpbar);
278  for (size_t k = 0; k < m_kk; k++) {
279  cpbar[k] *= GasConstant;
280  }
281 
282  // Check to see whether activity coefficients are temperature dependent. If
283  // they are, then calculate the their temperature derivatives and add them
284  // into the result.
285  double dAdT = dA_DebyedT_TP();
286  if (dAdT != 0.0) {
287  // Update the activity coefficients, This also update the internally
288  // stored molalities.
292  for (size_t k = 0; k < m_kk; k++) {
293  cpbar[k] -= (2.0 * RT() * m_dlnActCoeffMolaldT[k] +
295  }
296  }
297 }
298 
299 // -------------- Utilities -------------------------------
300 
301 //! Utility function to assign an integer value from a string for the
302 //! ElectrolyteSpeciesType field.
303 /*!
304  * @param estString input string that will be interpreted
305  */
306 static int interp_est(const std::string& estString)
307 {
308  if (caseInsensitiveEquals(estString, "solvent")) {
309  return cEST_solvent;
310  } else if (caseInsensitiveEquals(estString, "chargedspecies")) {
311  return cEST_chargedSpecies;
312  } else if (caseInsensitiveEquals(estString, "weakacidassociated")) {
313  return cEST_weakAcidAssociated;
314  } else if (caseInsensitiveEquals(estString, "strongacidassociated")) {
315  return cEST_strongAcidAssociated;
316  } else if (caseInsensitiveEquals(estString, "polarneutral")) {
317  return cEST_polarNeutral;
318  } else if (caseInsensitiveEquals(estString, "nonpolarneutral")) {
319  return cEST_nonpolarNeutral;
320  } else {
321  throw CanteraError("interp_est (DebyeHuckel)",
322  "Invalid electrolyte species type '{}'", estString);
323  }
324 }
325 
326 void DebyeHuckel::setDebyeHuckelModel(const std::string& model) {
327  if (model == "" || caseInsensitiveEquals(model, "Dilute_limit")) {
328  m_formDH = DHFORM_DILUTE_LIMIT;
329  } else if (caseInsensitiveEquals(model, "Bdot_with_variable_a")) {
330  m_formDH = DHFORM_BDOT_AK;
331  } else if (caseInsensitiveEquals(model, "Bdot_with_common_a")) {
332  m_formDH = DHFORM_BDOT_ACOMMON;
333  } else if (caseInsensitiveEquals(model, "Beta_ij")) {
334  m_formDH = DHFORM_BETAIJ;
335  m_Beta_ij.resize(m_kk, m_kk, 0.0);
336  } else if (caseInsensitiveEquals(model, "Pitzer_with_Beta_ij")) {
337  m_formDH = DHFORM_PITZER_BETAIJ;
338  m_Beta_ij.resize(m_kk, m_kk, 0.0);
339  } else {
340  throw CanteraError("DebyeHuckel::setDebyeHuckelModel",
341  "Unknown model '{}'", model);
342  }
343 }
344 
346 {
347  if (A < 0) {
348  m_form_A_Debye = A_DEBYE_WATER;
349  } else {
350  m_form_A_Debye = A_DEBYE_CONST;
351  m_A_Debye = A;
352  }
353 }
354 
355 void DebyeHuckel::setB_dot(double bdot)
356 {
357  if (m_formDH == DHFORM_BETAIJ || m_formDH == DHFORM_DILUTE_LIMIT ||
358  m_formDH == DHFORM_PITZER_BETAIJ) {
359  throw CanteraError("DebyeHuckel::setB_dot",
360  "B_dot entry in the wrong DH form");
361  }
362  // Set B_dot parameters for charged species
363  for (size_t k = 0; k < nSpecies(); k++) {
364  if (fabs(charge(k)) > 0.0001) {
365  m_B_Dot[k] = bdot;
366  } else {
367  m_B_Dot[k] = 0.0;
368  }
369  }
370 }
371 
373 {
374  for (size_t k = 0; k < m_kk; k++) {
375  if (std::isnan(m_Aionic[k])) {
376  m_Aionic[k] = value;
377  }
378  }
379 }
380 
381 void DebyeHuckel::setBeta(const std::string& sp1, const std::string& sp2,
382  double value)
383 {
384  size_t k1 = speciesIndex(sp1);
385  if (k1 == npos) {
386  throw CanteraError("DebyeHuckel::setBeta", "Species '{}' not found", sp1);
387  }
388  size_t k2 = speciesIndex(sp2);
389  if (k2 == npos) {
390  throw CanteraError("DebyeHuckel::setBeta", "Species '{}' not found", sp2);
391  }
392  m_Beta_ij(k1, k2) = value;
393  m_Beta_ij(k2, k1) = value;
394 }
395 
396 void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
397 {
398  if (id_.size() > 0) {
399  std::string idp = phaseNode.id();
400  if (idp != id_) {
401  throw CanteraError("DebyeHuckel::initThermoXML",
402  "phasenode and Id are incompatible");
403  }
404  }
405 
406  // Find the Thermo XML node
407  if (!phaseNode.hasChild("thermo")) {
408  throw CanteraError("DebyeHuckel::initThermoXML",
409  "no thermo XML node");
410  }
411  XML_Node& thermoNode = phaseNode.child("thermo");
412 
413  // Determine the form of the Debye-Huckel model, m_formDH. We will use
414  // this information to size arrays below. If there is no XML node named
415  // "activityCoefficients", assume that we are doing the extreme dilute
416  // limit assumption
417  if (thermoNode.hasChild("activityCoefficients")) {
418  setDebyeHuckelModel(thermoNode.child("activityCoefficients")["model"]);
419  } else {
420  setDebyeHuckelModel("Dilute_limit");
421  }
422 
423  // Go get all of the coefficients and factors in the activityCoefficients
424  // XML block
425  XML_Node* acNodePtr = 0;
426  if (thermoNode.hasChild("activityCoefficients")) {
427  XML_Node& acNode = thermoNode.child("activityCoefficients");
428  acNodePtr = &acNode;
429 
430  // Look for parameters for A_Debye
431  if (acNode.hasChild("A_Debye")) {
432  XML_Node* ss = acNode.findByName("A_Debye");
433  string modelString = ss->attrib("model");
434  if (modelString != "") {
435  if (caseInsensitiveEquals(modelString, "water")) {
436  setA_Debye(-1);
437  } else {
438  throw CanteraError("DebyeHuckel::initThermoXML",
439  "A_Debye Model \"" + modelString +
440  "\" is not known");
441  }
442  } else {
443  setA_Debye(getFloat(acNode, "A_Debye"));
444  }
445  }
446 
447  // Look for parameters for B_Debye
448  if (acNode.hasChild("B_Debye")) {
449  setB_Debye(getFloat(acNode, "B_Debye"));
450  }
451 
452  // Look for parameters for B_dot
453  if (acNode.hasChild("B_dot")) {
454  setB_dot(getFloat(acNode, "B_dot"));
455  }
456 
457  // Look for Parameters for the Maximum Ionic Strength
458  if (acNode.hasChild("maxIonicStrength")) {
459  setMaxIonicStrength(getFloat(acNode, "maxIonicStrength"));
460  }
461 
462  // Look for Helgeson Parameters
463  useHelgesonFixedForm(acNode.hasChild("UseHelgesonFixedForm"));
464 
465  // Look for parameters for the Ionic radius
466  if (acNode.hasChild("ionicRadius")) {
467  XML_Node& irNode = acNode.child("ionicRadius");
468 
469  double Afactor = 1.0;
470  if (irNode.hasAttrib("units")) {
471  std::string Aunits = irNode.attrib("units");
472  Afactor = toSI(Aunits);
473  }
474 
475  if (irNode.hasAttrib("default")) {
476  setDefaultIonicRadius(Afactor * fpValue(irNode.attrib("default")));
477  }
478 
479  // If the Debye-Huckel form is BDOT_AK, we can have separate values
480  // for the denominator's ionic size. -> That's how the activity
481  // coefficient is parameterized. In this case only do we allow the
482  // code to read in these parameters.
483  if (m_formDH == DHFORM_BDOT_AK) {
484  // Define a string-string map, and interpret the value of the
485  // XML element as binary pairs separated by colons, e.g.:
486  // Na+:3.0
487  // Cl-:4.0
488  // H+:9.0
489  // OH-:3.5
490  // Read them into the map.
491  map<string, string> m;
492  getMap(irNode, m);
493 
494  // Iterate over the map pairs, interpreting the first string as
495  // a species in the current phase. If no match is made, silently
496  // ignore the lack of agreement (HKM -> may be changed in the
497  // future).
498  for (const auto& b : m) {
499  size_t k = speciesIndex(b.first);
500  if (k != npos) {
501  m_Aionic[k] = fpValue(b.second) * Afactor;
502  }
503  }
504  }
505  }
506 
507  // Get the matrix of coefficients for the Beta binary interaction
508  // parameters. We assume here that this matrix is symmetric, so that we
509  // only have to input 1/2 of the values.
510  if (acNode.hasChild("DHBetaMatrix")) {
511  if (m_formDH == DHFORM_BETAIJ ||
512  m_formDH == DHFORM_PITZER_BETAIJ) {
513  XML_Node& irNode = acNode.child("DHBetaMatrix");
514  const vector<string>& sn = speciesNames();
515  getMatrixValues(irNode, sn, sn, m_Beta_ij, true, true);
516  } else {
517  throw CanteraError("DebyeHuckel::initThermoXML:",
518  "DHBetaMatrix found for wrong type");
519  }
520  }
521 
522  // Override stoichiometric Ionic Strength based on the phase definition
523  if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) {
524  XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
525  map<std::string, std::string> msIs;
526  getMap(sIsNode, msIs);
527  for (const auto& b : msIs) {
528  size_t kk = speciesIndex(b.first);
529  double val = fpValue(b.second);
530  m_speciesCharge_Stoich[kk] = val;
531  }
532  }
533  }
534 
535  // Override electrolyte species type based on the phase definition
536  if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) {
537  XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
538  map<std::string, std::string> msEST;
539  getMap(ESTNode, msEST);
540  for (const auto& b : msEST) {
541  size_t kk = speciesIndex(b.first);
542  std::string est = b.second;
543  if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
544  throw CanteraError("DebyeHuckel:initThermoXML",
545  "Bad electrolyte type: " + est);
546  }
547  }
548  }
549 
550  // Lastly set the state
551  if (phaseNode.hasChild("state")) {
552  XML_Node& stateNode = phaseNode.child("state");
553  setStateFromXML(stateNode);
554  }
555 }
556 
558 {
560  // Solvent
561  m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
562  if (m_waterSS) {
563  // Initialize the water property calculator. It will share the internal
564  // eos water calculator.
565  if (m_form_A_Debye == A_DEBYE_WATER) {
566  m_waterProps.reset(new WaterProps(m_waterSS));
567  }
568  } else if (dynamic_cast<PDSS_ConstVol*>(providePDSS(0)) == 0) {
569  throw CanteraError("DebyeHuckel::initThermo", "Solvent standard state"
570  " model must be WaterIAPWS or constant_incompressible.");
571  }
572 
573  // Solutes
574  for (size_t k = 1; k < nSpecies(); k++) {
575  if (dynamic_cast<PDSS_ConstVol*>(providePDSS(k)) == 0) {
576  throw CanteraError("DebyeHuckel::initThermo", "Solute standard"
577  " state model must be constant_incompressible.");
578  }
579  }
580 }
581 
582 double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const
583 {
584  double T = temperature();
585  double A;
586  if (tempArg != -1.0) {
587  T = tempArg;
588  }
589  double P = pressure();
590  if (presArg != -1.0) {
591  P = presArg;
592  }
593 
594  switch (m_form_A_Debye) {
595  case A_DEBYE_CONST:
596  A = m_A_Debye;
597  break;
598  case A_DEBYE_WATER:
599  A = m_waterProps->ADebye(T, P, 0);
600  m_A_Debye = A;
601  break;
602  default:
603  throw CanteraError("DebyeHuckel::A_Debye_TP", "shouldn't be here");
604  }
605  return A;
606 }
607 
608 double DebyeHuckel::dA_DebyedT_TP(double tempArg, double presArg) const
609 {
610  double T = temperature();
611  if (tempArg != -1.0) {
612  T = tempArg;
613  }
614  double P = pressure();
615  if (presArg != -1.0) {
616  P = presArg;
617  }
618  double dAdT;
619  switch (m_form_A_Debye) {
620  case A_DEBYE_CONST:
621  dAdT = 0.0;
622  break;
623  case A_DEBYE_WATER:
624  dAdT = m_waterProps->ADebye(T, P, 1);
625  break;
626  default:
627  throw CanteraError("DebyeHuckel::dA_DebyedT_TP", "shouldn't be here");
628  }
629  return dAdT;
630 }
631 
632 double DebyeHuckel::d2A_DebyedT2_TP(double tempArg, double presArg) const
633 {
634  double T = temperature();
635  if (tempArg != -1.0) {
636  T = tempArg;
637  }
638  double P = pressure();
639  if (presArg != -1.0) {
640  P = presArg;
641  }
642  double d2AdT2;
643  switch (m_form_A_Debye) {
644  case A_DEBYE_CONST:
645  d2AdT2 = 0.0;
646  break;
647  case A_DEBYE_WATER:
648  d2AdT2 = m_waterProps->ADebye(T, P, 2);
649  break;
650  default:
651  throw CanteraError("DebyeHuckel::d2A_DebyedT2_TP", "shouldn't be here");
652  }
653  return d2AdT2;
654 }
655 
656 double DebyeHuckel::dA_DebyedP_TP(double tempArg, double presArg) const
657 {
658  double T = temperature();
659  if (tempArg != -1.0) {
660  T = tempArg;
661  }
662  double P = pressure();
663  if (presArg != -1.0) {
664  P = presArg;
665  }
666  double dAdP;
667  switch (m_form_A_Debye) {
668  case A_DEBYE_CONST:
669  dAdP = 0.0;
670  break;
671  case A_DEBYE_WATER:
672  dAdP = m_waterProps->ADebye(T, P, 3);
673  break;
674  default:
675  throw CanteraError("DebyeHuckel::dA_DebyedP_TP", "shouldn't be here");
676  }
677  return dAdP;
678 }
679 
680 // ---------- Other Property Functions
681 
682 double DebyeHuckel::AionicRadius(int k) const
683 {
684  return m_Aionic[k];
685 }
686 
687 // ------------ Private and Restricted Functions ------------------
688 
689 bool DebyeHuckel::addSpecies(shared_ptr<Species> spec)
690 {
691  bool added = MolalityVPSSTP::addSpecies(spec);
692  if (added) {
693  m_lnActCoeffMolal.push_back(0.0);
694  m_dlnActCoeffMolaldT.push_back(0.0);
695  m_d2lnActCoeffMolaldT2.push_back(0.0);
696  m_dlnActCoeffMolaldP.push_back(0.0);
697  m_B_Dot.push_back(0.0);
698  m_tmpV.push_back(0.0);
699 
700  if (spec->extra.hasKey("ionic_radius")) {
701  m_Aionic.push_back(spec->extra["ionic_radius"].asDouble());
702  } else {
703  m_Aionic.push_back(NAN); // NAN will be replaced with default value
704  }
705 
706  // Guess electrolyte species type based on charge properties
707  int est = cEST_nonpolarNeutral;
708  double stoichCharge = spec->charge;
709  if (fabs(spec->charge) > 0.0001) {
710  est = cEST_chargedSpecies;
711  }
712  if (spec->extra.hasKey("weak_acid_charge")) {
713  stoichCharge = spec->extra["weak_acid_charge"].asDouble();
714  if (fabs(stoichCharge - spec->charge) > 0.0001) {
715  est = cEST_weakAcidAssociated;
716  }
717  }
718  m_speciesCharge_Stoich.push_back(stoichCharge);
719 
720  if (m_electrolyteSpeciesType.size() == 0) {
721  est = cEST_solvent; // species 0 is the solvent
722  }
723 
724  // Apply override of the electrolyte species type
725  if (spec->extra.hasKey("electrolyte_species_type")) {
726  est = interp_est(spec->extra["electrolyte_species_type"].asString());
727  }
728  m_electrolyteSpeciesType.push_back(est);
729  }
730  return added;
731 }
732 
733 double DebyeHuckel::_nonpolarActCoeff(double IionicMolality)
734 {
735  // These are coefficients to describe the increase in activity coeff for
736  // non-polar molecules due to the electrolyte becoming stronger (the so-
737  // called salt-out effect)
738  const static double npActCoeff[] = {0.1127, -0.01049, 1.545E-3};
739  double I2 = IionicMolality * IionicMolality;
740  double l10actCoeff =
741  npActCoeff[0] * IionicMolality +
742  npActCoeff[1] * I2 +
743  npActCoeff[2] * I2 * IionicMolality;
744  return pow(10.0 , l10actCoeff);
745 }
746 
748 {
749  const double a0 = 1.454;
750  const double b0 = 0.02236;
751  const double c0 = 9.380E-3;
752  const double d0 = -5.362E-4;
753  double Is = m_IionicMolalityStoich;
754  if (Is <= 0.0) {
755  return 0.0;
756  }
757  double Is2 = Is * Is;
758  double bhat = 1.0 + a0 * sqrt(Is);
759  double func = bhat - 2.0 * log(bhat) - 1.0/bhat;
760  double v1 = m_A_Debye / (a0 * a0 * a0 * Is) * func;
761  double oc = 1.0 - v1 + b0 * Is / 2.0 + 2.0 * c0 * Is2 / 3.0
762  + 3.0 * d0 * Is2 * Is / 4.0;
763  return oc;
764 }
765 
767 {
768  // Update the internally stored vector of molalities
769  calcMolalities();
770  double oc = _osmoticCoeffHelgesonFixedForm();
771  double sum = 0.0;
772  for (size_t k = 1; k < m_kk; k++) {
773  sum += std::max(m_molalities[k], 0.0);
774  }
775  if (sum > 2.0 * m_maxIionicStrength) {
776  sum = 2.0 * m_maxIionicStrength;
777  };
778  return - m_Mnaught * sum * oc;
779 }
780 
782 {
783  double z_k, zs_k1, zs_k2;
784 
785  // Update the internally stored vector of molalities
786  calcMolalities();
787 
788  // Calculate the apparent (real) ionic strength.
789  //
790  // Note this is not the stoichiometric ionic strengh, where reactions of
791  // ions forming neutral salts are ignorred in calculating the ionic
792  // strength.
793  m_IionicMolality = 0.0;
794  for (size_t k = 0; k < m_kk; k++) {
795  z_k = m_speciesCharge[k];
796  m_IionicMolality += m_molalities[k] * z_k * z_k;
797  }
798  m_IionicMolality /= 2.0;
800 
801  // Calculate the stoichiometric ionic charge
803  for (size_t k = 0; k < m_kk; k++) {
804  z_k = m_speciesCharge[k];
805  zs_k1 = m_speciesCharge_Stoich[k];
806  if (z_k == zs_k1) {
807  m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
808  } else {
809  zs_k2 = z_k - zs_k1;
810  m_IionicMolalityStoich += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
811  }
812  }
813  m_IionicMolalityStoich /= 2.0;
815 
816  // Possibly update the stored value of the Debye-Huckel parameter A_Debye
817  // This parameter appears on the top of the activity coefficient expression.
818  // It depends on T (and P), as it depends explicitly on the temperature.
819  // Also, the dielectric constant is usually a fairly strong function of T,
820  // also.
821  m_A_Debye = A_Debye_TP();
822 
823  // Calculate a safe value for the mole fraction of the solvent
824  double xmolSolvent = moleFraction(0);
825  xmolSolvent = std::max(8.689E-3, xmolSolvent);
826 
827  int est;
828  double ac_nonPolar = 1.0;
829  double numTmp = m_A_Debye * sqrt(m_IionicMolality);
830  double denomTmp = m_B_Debye * sqrt(m_IionicMolality);
831  double coeff;
832  double lnActivitySolvent = 0.0;
833  double tmp;
834  double tmpLn;
835  double y, yp1, sigma;
836  switch (m_formDH) {
837  case DHFORM_DILUTE_LIMIT:
838  for (size_t k = 0; k < m_kk; k++) {
839  z_k = m_speciesCharge[k];
840  m_lnActCoeffMolal[k] = - z_k * z_k * numTmp;
841  }
842  lnActivitySolvent =
843  (xmolSolvent - 1.0)/xmolSolvent +
844  2.0 / 3.0 * m_A_Debye * m_Mnaught *
846  break;
847 
848  case DHFORM_BDOT_AK:
849  ac_nonPolar = _nonpolarActCoeff(m_IionicMolality);
850  for (size_t k = 0; k < m_kk; k++) {
851  est = m_electrolyteSpeciesType[k];
852  if (est == cEST_nonpolarNeutral) {
853  m_lnActCoeffMolal[k] = log(ac_nonPolar);
854  } else {
855  z_k = m_speciesCharge[k];
856  m_lnActCoeffMolal[k] =
857  - z_k * z_k * numTmp / (1.0 + denomTmp * m_Aionic[k])
858  + log(10.0) * m_B_Dot[k] * m_IionicMolality;
859  }
860  }
861 
862  lnActivitySolvent = (xmolSolvent - 1.0)/xmolSolvent;
863  coeff = 2.0 / 3.0 * m_A_Debye * m_Mnaught
864  * sqrt(m_IionicMolality);
865  tmp = 0.0;
866  if (denomTmp > 0.0) {
867  for (size_t k = 0; k < m_kk; k++) {
868  if (k != 0 || m_Aionic[k] != 0.0) {
869  y = denomTmp * m_Aionic[k];
870  yp1 = y + 1.0;
871  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
872  z_k = m_speciesCharge[k];
873  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
874  }
875  }
876  }
877  lnActivitySolvent += coeff * tmp;
878  tmp = 0.0;
879  for (size_t k = 1; k < m_kk; k++) {
880  z_k = m_speciesCharge[k];
881  if (z_k != 0.0) {
882  tmp += m_B_Dot[k] * m_molalities[k];
883  }
884  }
885  lnActivitySolvent -=
886  m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
887 
888  // Special section to implement the Helgeson fixed form for the water
889  // brine activity coefficient.
891  lnActivitySolvent = _lnactivityWaterHelgesonFixedForm();
892  }
893  break;
894 
895  case DHFORM_BDOT_ACOMMON:
896  denomTmp *= m_Aionic[0];
897  for (size_t k = 0; k < m_kk; k++) {
898  z_k = m_speciesCharge[k];
899  m_lnActCoeffMolal[k] =
900  - z_k * z_k * numTmp / (1.0 + denomTmp)
901  + log(10.0) * m_B_Dot[k] * m_IionicMolality;
902  }
903  if (denomTmp > 0.0) {
904  y = denomTmp;
905  yp1 = y + 1.0;
906  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
907  } else {
908  sigma = 0.0;
909  }
910  lnActivitySolvent =
911  (xmolSolvent - 1.0)/xmolSolvent +
912  2.0 /3.0 * m_A_Debye * m_Mnaught *
913  m_IionicMolality * sqrt(m_IionicMolality) * sigma;
914  tmp = 0.0;
915  for (size_t k = 1; k < m_kk; k++) {
916  z_k = m_speciesCharge[k];
917  if (z_k != 0.0) {
918  tmp += m_B_Dot[k] * m_molalities[k];
919  }
920  }
921  lnActivitySolvent -=
922  m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
923  break;
924 
925  case DHFORM_BETAIJ:
926  denomTmp = m_B_Debye * m_Aionic[0];
927  denomTmp *= sqrt(m_IionicMolality);
928  lnActivitySolvent =
929  (xmolSolvent - 1.0)/xmolSolvent;
930 
931  for (size_t k = 1; k < m_kk; k++) {
932  z_k = m_speciesCharge[k];
933  m_lnActCoeffMolal[k] =
934  - z_k * z_k * numTmp / (1.0 + denomTmp);
935  for (size_t j = 0; j < m_kk; j++) {
936  double beta = m_Beta_ij.value(k, j);
937  m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
938  }
939  }
940  if (denomTmp > 0.0) {
941  y = denomTmp;
942  yp1 = y + 1.0;
943  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1));
944  } else {
945  sigma = 0.0;
946  }
947  lnActivitySolvent =
948  (xmolSolvent - 1.0)/xmolSolvent +
949  2.0 /3.0 * m_A_Debye * m_Mnaught *
950  m_IionicMolality * sqrt(m_IionicMolality) * sigma;
951  tmp = 0.0;
952  for (size_t k = 0; k < m_kk; k++) {
953  for (size_t j = 0; j < m_kk; j++) {
954  tmp +=
955  m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
956  }
957  }
958  lnActivitySolvent -= m_Mnaught * tmp;
959  break;
960 
961  case DHFORM_PITZER_BETAIJ:
962  denomTmp = m_B_Debye * sqrt(m_IionicMolality);
963  denomTmp *= m_Aionic[0];
964  numTmp = m_A_Debye * sqrt(m_IionicMolality);
965  tmpLn = log(1.0 + denomTmp);
966  for (size_t k = 1; k < m_kk; k++) {
967  z_k = m_speciesCharge[k];
968  m_lnActCoeffMolal[k] =
969  - z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
970  m_lnActCoeffMolal[k] +=
971  - 2.0 * z_k * z_k * m_A_Debye * tmpLn /
972  (3.0 * m_B_Debye * m_Aionic[0]);
973  for (size_t j = 0; j < m_kk; j++) {
974  m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
975  m_Beta_ij.value(k, j);
976  }
977  }
978  sigma = 1.0 / (1.0 + denomTmp);
979  lnActivitySolvent =
980  (xmolSolvent - 1.0)/xmolSolvent +
981  2.0 /3.0 * m_A_Debye * m_Mnaught *
982  m_IionicMolality * sqrt(m_IionicMolality) * sigma;
983  tmp = 0.0;
984  for (size_t k = 0; k < m_kk; k++) {
985  for (size_t j = 0; j < m_kk; j++) {
986  tmp +=
987  m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
988  }
989  }
990  lnActivitySolvent -= m_Mnaught * tmp;
991  break;
992 
993  default:
994  throw CanteraError("DebyeHuckel::s_update_lnMolalityActCoeff", "ERROR");
995  }
996 
997  // Above, we calculated the ln(activitySolvent). Translate that into the
998  // molar-based activity coefficient by dividing by the solvent mole
999  // fraction. Solvents are not on the molality scale.
1000  xmolSolvent = moleFraction(0);
1001  m_lnActCoeffMolal[0] = lnActivitySolvent - log(xmolSolvent);
1002 }
1003 
1005 {
1006  double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
1007  // First we store dAdT explicitly here
1008  double dAdT = dA_DebyedT_TP();
1009  if (dAdT == 0.0) {
1010  for (size_t k = 0; k < m_kk; k++) {
1011  m_dlnActCoeffMolaldT[k] = 0.0;
1012  }
1013  return;
1014  }
1015 
1016  // Calculate a safe value for the mole fraction of the solvent
1017  double xmolSolvent = moleFraction(0);
1018  xmolSolvent = std::max(8.689E-3, xmolSolvent);
1019  double sqrtI = sqrt(m_IionicMolality);
1020  double numdAdTTmp = dAdT * sqrtI;
1021  double denomTmp = m_B_Debye * sqrtI;
1022  double d_lnActivitySolvent_dT = 0;
1023 
1024  switch (m_formDH) {
1025  case DHFORM_DILUTE_LIMIT:
1026  for (size_t k = 1; k < m_kk; k++) {
1028  m_lnActCoeffMolal[k] * dAdT / m_A_Debye;
1029  }
1030  d_lnActivitySolvent_dT = 2.0 / 3.0 * dAdT * m_Mnaught *
1032  m_dlnActCoeffMolaldT[0] = d_lnActivitySolvent_dT;
1033  break;
1034 
1035  case DHFORM_BDOT_AK:
1036  for (size_t k = 0; k < m_kk; k++) {
1037  z_k = m_speciesCharge[k];
1039  - z_k * z_k * numdAdTTmp / (1.0 + denomTmp * m_Aionic[k]);
1040  }
1041 
1042  m_dlnActCoeffMolaldT[0] = 0.0;
1043  coeff = 2.0 / 3.0 * dAdT * m_Mnaught * sqrtI;
1044  tmp = 0.0;
1045  if (denomTmp > 0.0) {
1046  for (size_t k = 0; k < m_kk; k++) {
1047  y = denomTmp * m_Aionic[k];
1048  yp1 = y + 1.0;
1049  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1050  z_k = m_speciesCharge[k];
1051  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
1052  }
1053  }
1054  m_dlnActCoeffMolaldT[0] += coeff * tmp;
1055  break;
1056 
1057  case DHFORM_BDOT_ACOMMON:
1058  denomTmp *= m_Aionic[0];
1059  for (size_t k = 0; k < m_kk; k++) {
1060  z_k = m_speciesCharge[k];
1062  - z_k * z_k * numdAdTTmp / (1.0 + denomTmp);
1063  }
1064  if (denomTmp > 0.0) {
1065  y = denomTmp;
1066  yp1 = y + 1.0;
1067  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1068  } else {
1069  sigma = 0.0;
1070  }
1071  m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
1072  m_IionicMolality * sqrtI * sigma;
1073  break;
1074 
1075  case DHFORM_BETAIJ:
1076  denomTmp *= m_Aionic[0];
1077  for (size_t k = 1; k < m_kk; k++) {
1078  z_k = m_speciesCharge[k];
1079  m_dlnActCoeffMolaldT[k] = -z_k*z_k * numdAdTTmp / (1.0 + denomTmp);
1080  }
1081  if (denomTmp > 0.0) {
1082  y = denomTmp;
1083  yp1 = y + 1.0;
1084  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1085  } else {
1086  sigma = 0.0;
1087  }
1088  m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
1089  m_IionicMolality * sqrtI * sigma;
1090  break;
1091 
1092  case DHFORM_PITZER_BETAIJ:
1093  denomTmp *= m_Aionic[0];
1094  tmpLn = log(1.0 + denomTmp);
1095  for (size_t k = 1; k < m_kk; k++) {
1096  z_k = m_speciesCharge[k];
1098  - z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
1099  - 2.0 * z_k * z_k * dAdT * tmpLn / (m_B_Debye * m_Aionic[0]);
1100  m_dlnActCoeffMolaldT[k] /= 3.0;
1101  }
1102 
1103  sigma = 1.0 / (1.0 + denomTmp);
1104  m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
1105  m_IionicMolality * sqrtI * sigma;
1106  break;
1107 
1108  default:
1109  throw CanteraError("DebyeHuckel::s_update_dlnMolalityActCoeff_dT",
1110  "ERROR");
1111  }
1112 }
1113 
1115 {
1116  double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
1117  double dAdT = dA_DebyedT_TP();
1118  double d2AdT2 = d2A_DebyedT2_TP();
1119  if (d2AdT2 == 0.0 && dAdT == 0.0) {
1120  for (size_t k = 0; k < m_kk; k++) {
1121  m_d2lnActCoeffMolaldT2[k] = 0.0;
1122  }
1123  return;
1124  }
1125 
1126  // Calculate a safe value for the mole fraction of the solvent
1127  double xmolSolvent = moleFraction(0);
1128  xmolSolvent = std::max(8.689E-3, xmolSolvent);
1129  double sqrtI = sqrt(m_IionicMolality);
1130  double numd2AdT2Tmp = d2AdT2 * sqrtI;
1131  double denomTmp = m_B_Debye * sqrtI;
1132 
1133  switch (m_formDH) {
1134  case DHFORM_DILUTE_LIMIT:
1135  for (size_t k = 0; k < m_kk; k++) {
1137  m_lnActCoeffMolal[k] * d2AdT2 / m_A_Debye;
1138  }
1139  break;
1140 
1141  case DHFORM_BDOT_AK:
1142  for (size_t k = 0; k < m_kk; k++) {
1143  z_k = m_speciesCharge[k];
1145  - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp * m_Aionic[k]);
1146  }
1147 
1148  m_d2lnActCoeffMolaldT2[0] = 0.0;
1149  coeff = 2.0 / 3.0 * d2AdT2 * m_Mnaught * sqrtI;
1150  tmp = 0.0;
1151  if (denomTmp > 0.0) {
1152  for (size_t k = 0; k < m_kk; k++) {
1153  y = denomTmp * m_Aionic[k];
1154  yp1 = y + 1.0;
1155  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1156  z_k = m_speciesCharge[k];
1157  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
1158  }
1159  }
1160  m_d2lnActCoeffMolaldT2[0] += coeff * tmp;
1161  break;
1162 
1163  case DHFORM_BDOT_ACOMMON:
1164  denomTmp *= m_Aionic[0];
1165  for (size_t k = 0; k < m_kk; k++) {
1166  z_k = m_speciesCharge[k];
1168  - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp);
1169  }
1170  if (denomTmp > 0.0) {
1171  y = denomTmp;
1172  yp1 = y + 1.0;
1173  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1174  } else {
1175  sigma = 0.0;
1176  }
1177  m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
1178  m_IionicMolality * sqrtI * sigma;
1179  break;
1180 
1181  case DHFORM_BETAIJ:
1182  denomTmp *= m_Aionic[0];
1183  for (size_t k = 1; k < m_kk; k++) {
1184  z_k = m_speciesCharge[k];
1185  m_d2lnActCoeffMolaldT2[k] = -z_k*z_k * numd2AdT2Tmp / (1.0 + denomTmp);
1186  }
1187  if (denomTmp > 0.0) {
1188  y = denomTmp;
1189  yp1 = y + 1.0;
1190  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1));
1191  } else {
1192  sigma = 0.0;
1193  }
1194  m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
1195  m_IionicMolality * sqrtI * sigma;
1196  break;
1197 
1198  case DHFORM_PITZER_BETAIJ:
1199  denomTmp *= m_Aionic[0];
1200  tmpLn = log(1.0 + denomTmp);
1201  for (size_t k = 1; k < m_kk; k++) {
1202  z_k = m_speciesCharge[k];
1204  - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
1205  - 2.0 * z_k * z_k * d2AdT2 * tmpLn / (m_B_Debye * m_Aionic[0]);
1206  m_d2lnActCoeffMolaldT2[k] /= 3.0;
1207  }
1208 
1209  sigma = 1.0 / (1.0 + denomTmp);
1210  m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
1211  m_IionicMolality * sqrtI * sigma;
1212  break;
1213 
1214  default:
1215  throw CanteraError("DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2",
1216  "ERROR");
1217  }
1218 }
1219 
1221 {
1222  double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
1223  int est;
1224  double dAdP = dA_DebyedP_TP();
1225  if (dAdP == 0.0) {
1226  for (size_t k = 0; k < m_kk; k++) {
1227  m_dlnActCoeffMolaldP[k] = 0.0;
1228  }
1229  return;
1230  }
1231 
1232  // Calculate a safe value for the mole fraction of the solvent
1233  double xmolSolvent = moleFraction(0);
1234  xmolSolvent = std::max(8.689E-3, xmolSolvent);
1235  double sqrtI = sqrt(m_IionicMolality);
1236  double numdAdPTmp = dAdP * sqrtI;
1237  double denomTmp = m_B_Debye * sqrtI;
1238 
1239  switch (m_formDH) {
1240  case DHFORM_DILUTE_LIMIT:
1241  for (size_t k = 0; k < m_kk; k++) {
1243  m_lnActCoeffMolal[k] * dAdP / m_A_Debye;
1244  }
1245  break;
1246 
1247  case DHFORM_BDOT_AK:
1248  for (size_t k = 0; k < m_kk; k++) {
1249  est = m_electrolyteSpeciesType[k];
1250  if (est == cEST_nonpolarNeutral) {
1251  m_lnActCoeffMolal[k] = 0.0;
1252  } else {
1253  z_k = m_speciesCharge[k];
1255  - z_k * z_k * numdAdPTmp / (1.0 + denomTmp * m_Aionic[k]);
1256  }
1257  }
1258 
1259  m_dlnActCoeffMolaldP[0] = 0.0;
1260  coeff = 2.0 / 3.0 * dAdP * m_Mnaught * sqrtI;
1261  tmp = 0.0;
1262  if (denomTmp > 0.0) {
1263  for (size_t k = 0; k < m_kk; k++) {
1264  y = denomTmp * m_Aionic[k];
1265  yp1 = y + 1.0;
1266  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1267  z_k = m_speciesCharge[k];
1268  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
1269  }
1270  }
1271  m_dlnActCoeffMolaldP[0] += coeff * tmp;
1272  break;
1273 
1274  case DHFORM_BDOT_ACOMMON:
1275  denomTmp *= m_Aionic[0];
1276  for (size_t k = 0; k < m_kk; k++) {
1277  z_k = m_speciesCharge[k];
1279  - z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
1280  }
1281  if (denomTmp > 0.0) {
1282  y = denomTmp;
1283  yp1 = y + 1.0;
1284  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1285  } else {
1286  sigma = 0.0;
1287  }
1289  2.0 /3.0 * dAdP * m_Mnaught *
1290  m_IionicMolality * sqrtI * sigma;
1291  break;
1292 
1293  case DHFORM_BETAIJ:
1294  denomTmp *= m_Aionic[0];
1295  for (size_t k = 1; k < m_kk; k++) {
1296  z_k = m_speciesCharge[k];
1297  m_dlnActCoeffMolaldP[k] = - z_k*z_k * numdAdPTmp / (1.0 + denomTmp);
1298  }
1299  if (denomTmp > 0.0) {
1300  y = denomTmp;
1301  yp1 = y + 1.0;
1302  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
1303  } else {
1304  sigma = 0.0;
1305  }
1306  m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
1307  m_IionicMolality * sqrtI * sigma;
1308  break;
1309 
1310  case DHFORM_PITZER_BETAIJ:
1311  denomTmp *= m_Aionic[0];
1312  tmpLn = log(1.0 + denomTmp);
1313  for (size_t k = 1; k < m_kk; k++) {
1314  z_k = m_speciesCharge[k];
1316  - z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
1317  - 2.0 * z_k * z_k * dAdP * tmpLn
1318  / (m_B_Debye * m_Aionic[0]);
1319  m_dlnActCoeffMolaldP[k] /= 3.0;
1320  }
1321 
1322  sigma = 1.0 / (1.0 + denomTmp);
1323  m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
1324  m_IionicMolality * sqrtI * sigma;
1325  break;
1326 
1327  default:
1328  throw CanteraError("DebyeHuckel::s_update_dlnMolalityActCoeff_dP",
1329  "ERROR");
1330  }
1331 }
1332 
1333 }
void getMap(const XML_Node &node, std::map< std::string, std::string > &m)
This routine is used to interpret the value portions of XML elements that contain colon separated pai...
Definition: ctml.cpp:330
const int cEST_solvent
Electrolyte species type.
Definition: electrolytes.h:18
double m_IionicMolalityStoich
Stoichiometric ionic strength on the molality scale.
Definition: DebyeHuckel.h:1016
doublereal fpValue(const std::string &val)
Translate a string into one doublereal value.
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
Definition: DebyeHuckel.h:1035
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
Definition: DebyeHuckel.cpp:98
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
static double _nonpolarActCoeff(double IionicMolality)
Static function that implements the non-polar species salt-out modifications.
static int interp_est(const std::string &estString)
Utility function to assign an integer value from a string for the ElectrolyteSpeciesType field...
const XML_Node * findByName(const std::string &nm, int depth=100000) const
This routine carries out a recursive search for an XML node based on the name of the node...
Definition: xml.cpp:695
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
vector_fp m_dlnActCoeffMolaldT
Derivative of log act coeff wrt T.
Definition: DebyeHuckel.h:1133
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with &#39;v&#39;.
Definition: Array.h:112
size_t speciesIndex(const std::string &name) const
Returns the index of a species named &#39;name&#39; within the Phase object.
Definition: Phase.cpp:175
doublereal m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
void setDefaultIonicRadius(double value)
Set the default ionic radius [m] for each species.
doublereal toSI(const std::string &unit)
Return the conversion factor to convert unit std::string &#39;unit&#39; to SI units.
Definition: global.cpp:135
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:471
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
std::unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition: DebyeHuckel.h:1098
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
vector_fp m_speciesCharge
Vector of species charges. length m_kk.
Definition: Phase.h:799
Header file for a common definitions used in electrolytes thermodynamics.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
Definition: DebyeHuckel.cpp:80
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
virtual void initThermo()
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the species standard states at the cu...
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
STL namespace.
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:607
bool m_useHelgesonFixedForm
If true, then the fixed for of Helgeson&#39;s activity for water is used instead of the rigorous form obt...
Definition: DebyeHuckel.h:1013
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...
void s_update_lnMolalityActCoeff() const
Calculate the log activity coefficients.
void setDebyeHuckelModel(const std::string &form)
Set the DebyeHuckel parameterization form.
virtual bool addSpecies(shared_ptr< Species > spec)
void setBeta(const std::string &sp1, const std::string &sp2, double value)
Set the value for the beta interaction between species sp1 and sp2.
double _osmoticCoeffHelgesonFixedForm() const
Formula for the osmotic coefficient that occurs in the GWB.
bool hasAttrib(const std::string &a) const
Tests whether the current node has an attribute with a particular name.
Definition: xml.cpp:541
virtual void setDensity(const doublereal rho)
Set the internally stored density (gm/m^3) of the phase.
Headers for the DebyeHuckel ThermoPhase object, which models dilute electrolyte solutions (see Thermo...
const std::vector< std::string > & speciesNames() const
Return a const reference to the vector of species names.
Definition: Phase.cpp:197
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:614
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:748
DebyeHuckel()
Default Constructor.
Definition: DebyeHuckel.cpp:29
virtual doublereal molarVolume() const
Return the molar volume at standard state.
Definition: PDSS.cpp:72
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
Definition: DebyeHuckel.cpp:86
doublereal molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:590
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. ...
The WaterProps class is used to house several approximation routines for properties of water...
Definition: WaterProps.h:94
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:1123
vector_fp m_Aionic
a_k = Size of the ionic species in the DH formulation. units = meters
Definition: DebyeHuckel.h:999
Class for the liquid water pressure dependent standard state.
Definition: PDSS_Water.h:49
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition: DebyeHuckel.h:1006
int m_formDH
form of the Debye-Huckel parameterization used in the model.
Definition: DebyeHuckel.h:983
virtual void _updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
double m_A_Debye
Current value of the Debye Constant, A_Debye.
Definition: DebyeHuckel.h:1057
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
double AionicRadius(int k=0) const
Reports the ionic radius of the kth species.
vector_fp m_lnActCoeffMolal
Logarithm of the activity coefficients on the molality scale.
Definition: DebyeHuckel.h:1130
void s_update_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the activity coefficient.
bool caseInsensitiveEquals(const std::string &input, const std::string &test)
Case insensitive equality predicate.
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
PDSS_Water * m_waterSS
Pointer to the Water standard state object.
Definition: DebyeHuckel.h:1089
vector_fp m_B_Dot
Array of B_Dot values.
Definition: DebyeHuckel.h:1083
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
void s_update_d2lnMolalityActCoeff_dT2() const
Calculate the temperature 2nd derivative of the activity coefficient.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
virtual void setStateFromXML(const XML_Node &state)
Set equation of state parameter values from XML entries.
doublereal & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Definition: Array.h:231
void importPhase(XML_Node &phase, ThermoPhase *th)
Import a phase information into an empty ThermoPhase object.
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)) ...
void s_update_dlnMolalityActCoeff_dT() const
Calculation of temperature derivative of activity coefficient.
virtual void getActivities(doublereal *ac) const
Get the array of non-dimensional activities at the current solution temperature, pressure, and solution concentration.
virtual doublereal pressure() const
Returns the current pressure of the phase.
virtual void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:536
vector_int m_electrolyteSpeciesType
Vector containing the electrolyte species type.
Definition: DebyeHuckel.h:996
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized concentrations.
XML_Node & child(const size_t n) const
Return a changeable reference to the n&#39;th child of the current node.
Definition: xml.cpp:546
void setA_Debye(double A)
Set the A_Debye parameter.
vector_fp m_speciesCharge_Stoich
Stoichiometric species charge -> This is for calculations of the ionic strength which ignore ion-ion ...
Definition: DebyeHuckel.h:1115
void getMatrixValues(const XML_Node &node, const std::vector< std::string > &keyStringRow, const std::vector< std::string > &keyStringCol, Array2D &retnValues, const bool convert, const bool matrixSymmetric)
This function interprets the value portion of an XML element as a series of "Matrix ids and entries" ...
Definition: ctml.cpp:362
double m_IionicMolality
Current value of the ionic strength on the molality scale.
Definition: DebyeHuckel.h:1002
virtual void getMolalityActivityCoefficients(doublereal *acMolality) const
Get the array of non-dimensional molality-based activity coefficients at the current solution tempera...
vector_fp m_molalities
Current value of the molalities of the species in the phase.
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:126
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:500
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:661
vector_fp m_d2lnActCoeffMolaldT2
2nd Derivative of log act coeff wrt T
Definition: DebyeHuckel.h:1136
vector_fp m_tmpV
vector of size m_kk, used as a temporary holding area.
Definition: DebyeHuckel.h:1101
std::string id() const
Return the id attribute, if present.
Definition: xml.cpp:428
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
double m_B_Debye
Current value of the constant that appears in the denominator.
Definition: DebyeHuckel.h:1075
Contains declarations for string manipulation functions within Cantera.
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
Definition: ctml.cpp:164
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Definition: DebyeHuckel.cpp:92
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
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...
Declarations for the class PDSS_ConstVol (pressure dependent standard state) which handles calculatio...
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
size_t m_kk
Number of species in the phase.
Definition: Phase.h:788
vector_fp m_dlnActCoeffMolaldP
Derivative of log act coeff wrt P.
Definition: DebyeHuckel.h:1139
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
virtual doublereal density() const
Return the standard state density at standard state.
Definition: PDSS_Water.cpp:217
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
double _lnactivityWaterHelgesonFixedForm() const
Formula for the log of the water activity that occurs in the GWB.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:577
virtual bool addSpecies(shared_ptr< Species > spec)
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of the species standard states at the current T and P of the solution...
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
Definition: Phase.h:622
double m_densWaterSS
Storage for the density of water&#39;s standard state.
Definition: DebyeHuckel.h:1095
virtual void setMolarDensity(const doublereal conc)
Set the internally stored molar density (kmol/m^3) of the phase.