Cantera  2.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  * Copyright (2006) Sandia Corporation. Under the terms of
12  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
13  * U.S. Government retains certain rights in this software.
14  */
15 
16 #include "cantera/thermo/DebyeHuckel.h"
17 #include "cantera/thermo/ThermoFactory.h"
18 #include "cantera/thermo/WaterProps.h"
19 #include "cantera/thermo/PDSS_Water.h"
20 
21 #include "cantera/base/stringUtils.h"
22 
23 #include <cstring>
24 #include <cstdlib>
25 #include <fstream>
26 
27 using namespace std;
28 using namespace ctml;
29 
30 namespace Cantera
31 {
32 
33 /*
34  * Default constructor
35  */
36 DebyeHuckel::DebyeHuckel() :
37  MolalityVPSSTP(),
38  m_formDH(DHFORM_DILUTE_LIMIT),
39  m_formGC(2),
40  m_IionicMolality(0.0),
41  m_maxIionicStrength(30.0),
42  m_useHelgesonFixedForm(false),
43  m_IionicMolalityStoich(0.0),
44  m_form_A_Debye(A_DEBYE_CONST),
45  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
46  m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
47  m_waterSS(0),
48  m_densWaterSS(1000.),
49  m_waterProps(0)
50 {
51 
52  m_npActCoeff.resize(3);
53  m_npActCoeff[0] = 0.1127;
54  m_npActCoeff[1] = -0.01049;
55  m_npActCoeff[2] = 1.545E-3;
56 }
57 
58 /*
59  * Working constructors
60  *
61  * The two constructors below are the normal way
62  * the phase initializes itself. They are shells that call
63  * the routine initThermo(), with a reference to the
64  * XML database to get the info for the phase.
65  */
66 DebyeHuckel::DebyeHuckel(std::string inputFile, std::string id) :
67  MolalityVPSSTP(),
68  m_formDH(DHFORM_DILUTE_LIMIT),
69  m_formGC(2),
70  m_IionicMolality(0.0),
71  m_maxIionicStrength(30.0),
72  m_useHelgesonFixedForm(false),
73  m_IionicMolalityStoich(0.0),
74  m_form_A_Debye(A_DEBYE_CONST),
75  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
76  m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
77  m_waterSS(0),
78  m_densWaterSS(1000.),
79  m_waterProps(0)
80 {
81  m_npActCoeff.resize(3);
82  m_npActCoeff[0] = 0.1127;
83  m_npActCoeff[1] = -0.01049;
84  m_npActCoeff[2] = 1.545E-3;
85  constructPhaseFile(inputFile, id);
86 }
87 
88 DebyeHuckel::DebyeHuckel(XML_Node& phaseRoot, std::string id) :
89  MolalityVPSSTP(),
90  m_formDH(DHFORM_DILUTE_LIMIT),
91  m_formGC(2),
92  m_IionicMolality(0.0),
93  m_maxIionicStrength(3.0),
94  m_useHelgesonFixedForm(false),
95  m_IionicMolalityStoich(0.0),
96  m_form_A_Debye(A_DEBYE_CONST),
97  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
98  m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
99  m_waterSS(0),
100  m_densWaterSS(1000.),
101  m_waterProps(0)
102 {
103  m_npActCoeff.resize(3);
104  m_npActCoeff[0] = 0.1127;
105  m_npActCoeff[1] = -0.01049;
106  m_npActCoeff[2] = 1.545E-3;
107  constructPhaseXML(phaseRoot, id);
108 }
109 
110 /*
111  * Copy Constructor:
112  *
113  * Note this stuff will not work until the underlying phase
114  * has a working copy constructor
115  */
116 DebyeHuckel::DebyeHuckel(const DebyeHuckel& b) :
117  MolalityVPSSTP(),
118  m_formDH(DHFORM_DILUTE_LIMIT),
119  m_formGC(2),
120  m_IionicMolality(0.0),
121  m_maxIionicStrength(30.0),
122  m_useHelgesonFixedForm(false),
123  m_IionicMolalityStoich(0.0),
124  m_form_A_Debye(A_DEBYE_CONST),
125  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
126  m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
127  m_waterSS(0),
128  m_densWaterSS(1000.),
129  m_waterProps(0)
130 {
131  /*
132  * Use the assignment operator to do the brunt
133  * of the work for the copy constructor.
134  */
135  *this = b;
136 }
137 
138 /*
139  * operator=()
140  *
141  * Note this stuff will not work until the underlying phase
142  * has a working assignment operator
143  */
144 DebyeHuckel& DebyeHuckel::
145 operator=(const DebyeHuckel& b)
146 {
147  if (&b != this) {
148  MolalityVPSSTP::operator=(b);
149  m_formDH = b.m_formDH;
150  m_formGC = b.m_formGC;
151  m_Aionic = b.m_Aionic;
152  m_npActCoeff = b.m_npActCoeff;
153  m_IionicMolality = b.m_IionicMolality;
154  m_maxIionicStrength = b.m_maxIionicStrength;
155  m_useHelgesonFixedForm= b.m_useHelgesonFixedForm;
156  m_IionicMolalityStoich= b.m_IionicMolalityStoich;
157  m_form_A_Debye = b.m_form_A_Debye;
158  m_A_Debye = b.m_A_Debye;
159  m_B_Debye = b.m_B_Debye;
160  m_B_Dot = b.m_B_Dot;
161  m_npActCoeff = b.m_npActCoeff;
162 
163  // This is an internal shallow copy of the PDSS_Water pointer
164  m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0)) ;
165  if (!m_waterSS) {
166  throw CanteraError("DebyHuckel::operator=()", "Dynamic cast to waterPDSS failed");
167  }
168 
169  m_densWaterSS = b.m_densWaterSS;
170 
171  if (m_waterProps) {
172  delete m_waterProps;
173  m_waterProps = 0;
174  }
175  if (b.m_waterProps) {
176  m_waterProps = new WaterProps(m_waterSS);
177  }
178 
179  m_expg0_RT = b.m_expg0_RT;
180  m_pe = b.m_pe;
181  m_pp = b.m_pp;
182  m_tmpV = b.m_tmpV;
183  m_speciesCharge_Stoich= b.m_speciesCharge_Stoich;
184  m_Beta_ij = b.m_Beta_ij;
185  m_lnActCoeffMolal = b.m_lnActCoeffMolal;
186  m_d2lnActCoeffMolaldT2= b.m_d2lnActCoeffMolaldT2;
187  }
188  return *this;
189 }
190 
191 
192 /*
193  * ~DebyeHuckel(): (virtual)
194  *
195  * Destructor for DebyeHuckel. Release objects that
196  * it owns.
197  */
198 DebyeHuckel::~DebyeHuckel()
199 {
200  if (m_waterProps) {
201  delete m_waterProps;
202  m_waterProps = 0;
203  }
204 }
205 
206 /*
207  * duplMyselfAsThermoPhase():
208  *
209  * This routine operates at the ThermoPhase level to
210  * duplicate the current object. It uses the copy constructor
211  * defined above.
212  */
213 ThermoPhase* DebyeHuckel::duplMyselfAsThermoPhase() const
214 {
215  DebyeHuckel* mtp = new DebyeHuckel(*this);
216  return (ThermoPhase*) mtp;
217 }
218 
219 /*
220  * Equation of state type flag. The base class returns
221  * zero. Subclasses should define this to return a unique
222  * non-zero value. Constants defined for this purpose are
223  * listed in mix_defs.h.
224  */
225 int DebyeHuckel::eosType() const
226 {
227  int res;
228  switch (m_formGC) {
229  case 0:
230  res = cDebyeHuckel0;
231  break;
232  case 1:
233  res = cDebyeHuckel1;
234  break;
235  case 2:
236  res = cDebyeHuckel2;
237  break;
238  default:
239  throw CanteraError("eosType", "Unknown type");
240  break;
241  }
242  return res;
243 }
244 
245 //
246 // -------- Molar Thermodynamic Properties of the Solution ---------------
247 //
248 /*
249  * Molar enthalpy of the solution. Units: J/kmol.
250  */
251 doublereal DebyeHuckel::enthalpy_mole() const
252 {
253  getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
254  return mean_X(DATA_PTR(m_tmpV));
255 }
256 
257 /*
258  * Molar internal energy of the solution. Units: J/kmol.
259  *
260  * This is calculated from the soln enthalpy and then
261  * subtracting pV.
262  */
263 doublereal DebyeHuckel::intEnergy_mole() const
264 {
265  double hh = enthalpy_mole();
266  double pres = pressure();
267  double molarV = 1.0/molarDensity();
268  double uu = hh - pres * molarV;
269  return uu;
270 }
271 
272 /*
273  * Molar soln entropy at constant pressure. Units: J/kmol/K.
274  *
275  * This is calculated from the partial molar entropies.
276  */
277 doublereal DebyeHuckel::entropy_mole() const
278 {
279  getPartialMolarEntropies(DATA_PTR(m_tmpV));
280  return mean_X(DATA_PTR(m_tmpV));
281 }
282 
283 // Molar Gibbs function. Units: J/kmol.
284 doublereal DebyeHuckel::gibbs_mole() const
285 {
286  getChemPotentials(DATA_PTR(m_tmpV));
287  return mean_X(DATA_PTR(m_tmpV));
288 }
289 
290 /*
291  * Molar heat capacity at constant pressure. Units: J/kmol/K.
292  *
293  * Returns the solution heat capacition at constant pressure.
294  * This is calculated from the partial molar heat capacities.
295  */
296 doublereal DebyeHuckel::cp_mole() const
297 {
298  getPartialMolarCp(DATA_PTR(m_tmpV));
299  double val = mean_X(DATA_PTR(m_tmpV));
300  return val;
301 }
302 
303 /// Molar heat capacity at constant volume. Units: J/kmol/K.
304 doublereal DebyeHuckel::cv_mole() const
305 {
306  //getPartialMolarCv(m_tmpV.begin());
307  //return mean_X(m_tmpV.begin());
308  err("not implemented");
309  return 0.0;
310 }
311 
312 //
313 // ------- Mechanical Equation of State Properties ------------------------
314 //
315 
316 /*
317  * Pressure. Units: Pa.
318  * For this incompressible system, we return the internally stored
319  * independent value of the pressure.
320  */
321 doublereal DebyeHuckel::pressure() const
322 {
323  return m_Pcurrent;
324 }
325 
326 void DebyeHuckel::setPressure(doublereal p)
327 {
328  setState_TP(temperature(), p);
329 }
330 
331 void DebyeHuckel::setState_TP(doublereal t, doublereal p)
332 {
333 
334  Phase::setTemperature(t);
335  /*
336  * Store the current pressure
337  */
338  m_Pcurrent = p;
339 
340  /*
341  * update the standard state thermo
342  * -> This involves calling the water function and setting the pressure
343  */
344  _updateStandardStateThermo();
345 
346  /*
347  * Calculate all of the other standard volumes
348  * -> note these are constant for now
349  */
350  calcDensity();
351 }
352 
353 /*
354  * Calculate the density of the mixture using the partial
355  * molar volumes and mole fractions as input
356  *
357  * The formula for this is
358  *
359  * \f[
360  * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}}
361  * \f]
362  *
363  * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are
364  * the molecular weights, and \f$V_k\f$ are the pure species
365  * molar volumes.
366  *
367  * Note, the basis behind this formula is that in an ideal
368  * solution the partial molar volumes are equal to the pure
369  * species molar volumes. We have additionally specified
370  * in this class that the pure species molar volumes are
371  * independent of temperature and pressure.
372  *
373  */
374 void DebyeHuckel::calcDensity()
375 {
376  if (m_waterSS) {
377 
378  /*
379  * Store the internal density of the water SS.
380  * Note, we would have to do this for all other
381  * species if they had pressure dependent properties.
382  */
383  m_densWaterSS = m_waterSS->density();
384  }
385  double* vbar = &m_pp[0];
386  getPartialMolarVolumes(vbar);
387  double* x = &m_tmpV[0];
388  getMoleFractions(x);
389  doublereal vtotal = 0.0;
390  for (size_t i = 0; i < m_kk; i++) {
391  vtotal += vbar[i] * x[i];
392  }
393  doublereal dd = meanMolecularWeight() / vtotal;
394  Phase::setDensity(dd);
395 }
396 
397 
398 /*
399  * The isothermal compressibility. Units: 1/Pa.
400  * The isothermal compressibility is defined as
401  * \f[
402  * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
403  * \f]
404  *
405  * It's equal to zero for this model, since the molar volume
406  * doesn't change with pressure or temperature.
407  */
408 doublereal DebyeHuckel::isothermalCompressibility() const
409 {
410  throw CanteraError("DebyeHuckel::isothermalCompressibility",
411  "unimplemented");
412  return 0.0;
413 }
414 
415 /*
416  * The thermal expansion coefficient. Units: 1/K.
417  * The thermal expansion coefficient is defined as
418  *
419  * \f[
420  * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
421  * \f]
422  *
423  * It's equal to zero for this model, since the molar volume
424  * doesn't change with pressure or temperature.
425  */
426 doublereal DebyeHuckel::thermalExpansionCoeff() const
427 {
428  throw CanteraError("DebyeHuckel::thermalExpansionCoeff",
429  "unimplemented");
430  return 0.0;
431 }
432 
433 /*
434  * Overwritten setDensity() function is necessary because the
435  * density is not an independent variable.
436  *
437  * This function will now throw an error condition
438  *
439  * @internal May have to adjust the strategy here to make
440  * the eos for these materials slightly compressible, in order
441  * to create a condition where the density is a function of
442  * the pressure.
443  *
444  * This function will now throw an error condition.
445  *
446  * NOTE: This is an overwritten function from the State.h
447  * class
448  */
449 void DebyeHuckel::setDensity(doublereal rho)
450 {
451  double dens = density();
452  if (rho != dens) {
453  throw CanteraError("Idea;MolalSoln::setDensity",
454  "Density is not an independent variable");
455  }
456 }
457 
458 /*
459  * Overwritten setMolarDensity() function is necessary because the
460  * density is not an independent variable.
461  *
462  * This function will now throw an error condition.
463  *
464  * NOTE: This is a virtual function, overwritten function from the State.h
465  * class
466  */
467 void DebyeHuckel::setMolarDensity(const doublereal conc)
468 {
469  double concI = molarDensity();
470  if (conc != concI) {
471  throw CanteraError("Idea;MolalSoln::setMolarDensity",
472  "molarDensity/density is not an independent variable");
473  }
474 }
475 
476 /*
477  * Overwritten setTemperature(double) from State.h. This
478  * function sets the temperature, and makes sure that
479  * the value propagates to underlying objects.
480  */
481 void DebyeHuckel::setTemperature(const doublereal temp)
482 {
483  setState_TP(temp, m_Pcurrent);
484 }
485 
486 
487 //
488 // ------- Activities and Activity Concentrations
489 //
490 
491 /*
492  * This method returns an array of generalized concentrations
493  * \f$ C_k\f$ that are defined such that
494  * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
495  * is a standard concentration
496  * defined below. These generalized concentrations are used
497  * by kinetics manager classes to compute the forward and
498  * reverse rates of elementary reactions.
499  *
500  * @param c Array of generalized concentrations. The
501  * units depend upon the implementation of the
502  * reaction rate expressions within the phase.
503  */
504 void DebyeHuckel::getActivityConcentrations(doublereal* c) const
505 {
506  double c_solvent = standardConcentration();
507  getActivities(c);
508  for (size_t k = 0; k < m_kk; k++) {
509  c[k] *= c_solvent;
510  }
511 }
512 
513 /*
514  * The standard concentration \f$ C^0_k \f$ used to normalize
515  * the generalized concentration. In many cases, this quantity
516  * will be the same for all species in a phase - for example,
517  * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
518  * reason, this method returns a single value, instead of an
519  * array. However, for phases in which the standard
520  * concentration is species-specific (e.g. surface species of
521  * different sizes), this method may be called with an
522  * optional parameter indicating the species.
523  *
524  * For the time being we will use the concentration of pure
525  * solvent for the the standard concentration of all species.
526  * This has the effect of making reaction rates
527  * based on the molality of species proportional to the
528  * molality of the species.
529  */
530 doublereal DebyeHuckel::standardConcentration(size_t k) const
531 {
532  double mvSolvent = m_speciesSize[m_indexSolvent];
533  return 1.0 / mvSolvent;
534 }
535 
536 /*
537  * Returns the natural logarithm of the standard
538  * concentration of the kth species
539  */
540 doublereal DebyeHuckel::logStandardConc(size_t k) const
541 {
542  double c_solvent = standardConcentration(k);
543  return log(c_solvent);
544 }
545 
546 /*
547  * Returns the units of the standard and general concentrations
548  * Note they have the same units, as their divisor is
549  * defined to be equal to the activity of the kth species
550  * in the solution, which is unitless.
551  *
552  * This routine is used in print out applications where the
553  * units are needed. Usually, MKS units are assumed throughout
554  * the program and in the XML input files.
555  *
556  * On return uA contains the powers of the units (MKS assumed)
557  * of the standard concentrations and generalized concentrations
558  * for the kth species.
559  *
560  * uA[0] = kmol units - default = 1
561  * uA[1] = m units - default = -nDim(), the number of spatial
562  * dimensions in the Phase class.
563  * uA[2] = kg units - default = 0;
564  * uA[3] = Pa(pressure) units - default = 0;
565  * uA[4] = Temperature units - default = 0;
566  * uA[5] = time units - default = 0
567  */
568 void DebyeHuckel::getUnitsStandardConc(double* uA, int k, int sizeUA) const
569 {
570  for (int i = 0; i < sizeUA; i++) {
571  if (i == 0) {
572  uA[0] = 1.0;
573  }
574  if (i == 1) {
575  uA[1] = -int(nDim());
576  }
577  if (i == 2) {
578  uA[2] = 0.0;
579  }
580  if (i == 3) {
581  uA[3] = 0.0;
582  }
583  if (i == 4) {
584  uA[4] = 0.0;
585  }
586  if (i == 5) {
587  uA[5] = 0.0;
588  }
589  }
590 }
591 
592 
593 /*
594  * Get the array of non-dimensional activities at
595  * the current solution temperature, pressure, and
596  * solution concentration.
597  * (note solvent activity coefficient is on the molar scale).
598  *
599  */
600 void DebyeHuckel::getActivities(doublereal* ac) const
601 {
602  _updateStandardStateThermo();
603  /*
604  * Update the molality array, m_molalities()
605  * This requires an update due to mole fractions
606  */
607  s_update_lnMolalityActCoeff();
608  for (size_t k = 0; k < m_kk; k++) {
609  if (k != m_indexSolvent) {
610  ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
611  }
612  }
613  double xmolSolvent = moleFraction(m_indexSolvent);
614  ac[m_indexSolvent] =
615  exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent;
616 }
617 
618 /*
619  * getMolalityActivityCoefficients() (virtual, const)
620  *
621  * Get the array of non-dimensional Molality based
622  * activity coefficients at
623  * the current solution temperature, pressure, and
624  * solution concentration.
625  * (note solvent activity coefficient is on the molar scale).
626  *
627  * Note, most of the work is done in an internal private routine
628  */
629 void DebyeHuckel::
630 getMolalityActivityCoefficients(doublereal* acMolality) const
631 {
632  _updateStandardStateThermo();
633  A_Debye_TP(-1.0, -1.0);
634  s_update_lnMolalityActCoeff();
635  copy(m_lnActCoeffMolal.begin(), m_lnActCoeffMolal.end(), acMolality);
636  for (size_t k = 0; k < m_kk; k++) {
637  acMolality[k] = exp(acMolality[k]);
638  }
639 }
640 
641 //
642 // ------ Partial Molar Properties of the Solution -----------------
643 //
644 /*
645  * Get the species chemical potentials. Units: J/kmol.
646  *
647  * This function returns a vector of chemical potentials of the
648  * species in solution.
649  *
650  * \f[
651  * \mu_k = \mu^{o}_k(T,P) + R T ln(m_k)
652  * \f]
653  *
654  * \f[
655  * \mu_solvent = \mu^{o}_solvent(T,P) +
656  * R T ((X_solvent - 1.0) / X_solvent)
657  * \f]
658  */
659 void DebyeHuckel::getChemPotentials(doublereal* mu) const
660 {
661  double xx;
662  const double xxSmall = 1.0E-150;
663  /*
664  * First get the standard chemical potentials in
665  * molar form.
666  * -> this requires updates of standard state as a function
667  * of T and P
668  */
669  getStandardChemPotentials(mu);
670  /*
671  * Update the activity coefficients
672  * This also updates the internal molality array.
673  */
674  s_update_lnMolalityActCoeff();
675  doublereal RT = GasConstant * temperature();
676  double xmolSolvent = moleFraction(m_indexSolvent);
677  for (size_t k = 0; k < m_kk; k++) {
678  if (m_indexSolvent != k) {
679  xx = std::max(m_molalities[k], xxSmall);
680  mu[k] += RT * (log(xx) + m_lnActCoeffMolal[k]);
681  }
682  }
683  xx = std::max(xmolSolvent, xxSmall);
684  mu[m_indexSolvent] +=
685  RT * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]);
686 }
687 
688 
689 /*
690  * Returns an array of partial molar enthalpies for the species
691  * in the mixture.
692  * Units (J/kmol)
693  *
694  * We calculate this quantity partially from the relation and
695  * partially by calling the standard state enthalpy function.
696  *
697  * hbar_i = - T**2 * d(chemPot_i/T)/dT
698  *
699  * We calculate
700  */
701 void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const
702 {
703  /*
704  * Get the nondimensional standard state enthalpies
705  */
706  getEnthalpy_RT(hbar);
707  /*
708  * Dimensionalize it.
709  */
710  double T = temperature();
711  double RT = GasConstant * T;
712  for (size_t k = 0; k < m_kk; k++) {
713  hbar[k] *= RT;
714  }
715  /*
716  * Check to see whether activity coefficients are temperature
717  * dependent. If they are, then calculate the their temperature
718  * derivatives and add them into the result.
719  */
720  double dAdT = dA_DebyedT_TP();
721  if (dAdT != 0.0) {
722  /*
723  * Update the activity coefficients, This also update the
724  * internally stored molalities.
725  */
726  s_update_lnMolalityActCoeff();
727  s_update_dlnMolalityActCoeff_dT();
728  double RTT = GasConstant * T * T;
729  for (size_t k = 0; k < m_kk; k++) {
730  hbar[k] -= RTT * m_dlnActCoeffMolaldT[k];
731  }
732  }
733 }
734 
735 /*
736  *
737  * getPartialMolarEntropies() (virtual, const)
738  *
739  * Returns an array of partial molar entropies of the species in the
740  * solution. Units: J/kmol.
741  *
742  * Maxwell's equations provide an insight in how to calculate this
743  * (p.215 Smith and Van Ness)
744  *
745  * d(chemPot_i)/dT = -sbar_i
746  *
747  * For this phase, the partial molar entropies are equal to the
748  * SS species entropies plus the ideal solution contribution.following
749  * contribution:
750  * \f[
751  * \bar s_k(T,P) = \hat s^0_k(T) - R log(M0 * molality[k])
752  * \f]
753  * \f[
754  * \bar s_solvent(T,P) = \hat s^0_solvent(T)
755  * - R ((xmolSolvent - 1.0) / xmolSolvent)
756  * \f]
757  *
758  * The reference-state pure-species entropies,\f$ \hat s^0_k(T) \f$,
759  * at the reference pressure, \f$ P_{ref} \f$, are computed by the
760  * species thermodynamic
761  * property manager. They are polynomial functions of temperature.
762  * @see SpeciesThermo
763  */
764 void DebyeHuckel::
765 getPartialMolarEntropies(doublereal* sbar) const
766 {
767  /*
768  * Get the standard state entropies at the temperature
769  * and pressure of the solution.
770  */
771  getEntropy_R(sbar);
772  /*
773  * Dimensionalize the entropies
774  */
775  doublereal R = GasConstant;
776  for (size_t k = 0; k < m_kk; k++) {
777  sbar[k] *= R;
778  }
779  /*
780  * Update the activity coefficients, This also update the
781  * internally stored molalities.
782  */
783  s_update_lnMolalityActCoeff();
784  /*
785  * First we will add in the obvious dependence on the T
786  * term out front of the log activity term
787  */
788  doublereal mm;
789  for (size_t k = 0; k < m_kk; k++) {
790  if (k != m_indexSolvent) {
791  mm = std::max(SmallNumber, m_molalities[k]);
792  sbar[k] -= R * (log(mm) + m_lnActCoeffMolal[k]);
793  }
794  }
795  double xmolSolvent = moleFraction(m_indexSolvent);
796  mm = std::max(SmallNumber, xmolSolvent);
797  sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
798  /*
799  * Check to see whether activity coefficients are temperature
800  * dependent. If they are, then calculate the their temperature
801  * derivatives and add them into the result.
802  */
803  double dAdT = dA_DebyedT_TP();
804  if (dAdT != 0.0) {
805  s_update_dlnMolalityActCoeff_dT();
806  double RT = R * temperature();
807  for (size_t k = 0; k < m_kk; k++) {
808  sbar[k] -= RT * m_dlnActCoeffMolaldT[k];
809  }
810  }
811 }
812 
813 /*
814  * getPartialMolarVolumes() (virtual, const)
815  *
816  * returns an array of partial molar volumes of the species
817  * in the solution. Units: m^3 kmol-1.
818  *
819  * For this solution, the partial molar volumes are normally
820  * equal to theconstant species molar volumes, except
821  * when the activity coefficients depend on pressure.
822  *
823  * The general relation is
824  *
825  * vbar_i = d(chemPot_i)/dP at const T, n
826  *
827  * = V0_i + d(Gex)/dP)_T,M
828  *
829  * = V0_i + RT d(lnActCoeffi)dP _T,M
830  *
831  */
832 void DebyeHuckel::getPartialMolarVolumes(doublereal* vbar) const
833 {
834  getStandardVolumes(vbar);
835  /*
836  * Update the derivatives wrt the activity coefficients.
837  */
838  s_update_lnMolalityActCoeff();
839  s_update_dlnMolalityActCoeff_dP();
840  double T = temperature();
841  double RT = GasConstant * T;
842  for (size_t k = 0; k < m_kk; k++) {
843  vbar[k] += RT * m_dlnActCoeffMolaldP[k];
844  }
845 }
846 
847 /*
848  * Partial molar heat capacity of the solution:
849  * The kth partial molar heat capacity is equal to
850  * the temperature derivative of the partial molar
851  * enthalpy of the kth species in the solution at constant
852  * P and composition (p. 220 Smith and Van Ness).
853  *
854  * Cp = -T d2(chemPot_i)/dT2
855  */
856 void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const
857 {
858  /*
859  * Get the nondimensional gibbs standard state of the
860  * species at the T and P of the solution.
861  */
862  getCp_R(cpbar);
863 
864  for (size_t k = 0; k < m_kk; k++) {
865  cpbar[k] *= GasConstant;
866  }
867 
868  /*
869  * Check to see whether activity coefficients are temperature
870  * dependent. If they are, then calculate the their temperature
871  * derivatives and add them into the result.
872  */
873  double dAdT = dA_DebyedT_TP();
874  if (dAdT != 0.0) {
875  /*
876  * Update the activity coefficients, This also update the
877  * internally stored molalities.
878  */
879  s_update_lnMolalityActCoeff();
880  s_update_dlnMolalityActCoeff_dT();
881  s_update_d2lnMolalityActCoeff_dT2();
882  double T = temperature();
883  double RT = GasConstant * T;
884  double RTT = RT * T;
885  for (size_t k = 0; k < m_kk; k++) {
886  cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT[k] +
887  RTT * m_d2lnActCoeffMolaldT2[k]);
888  }
889  }
890 }
891 
892 
893 
894 
895 /*
896  * -------------- Utilities -------------------------------
897  */
898 
899 /*
900  * Initialization routine for a DebyeHuckel phase.
901  *
902  * This is a virtual routine. This routine will call initThermo()
903  * for the parent class as well.
904  */
905 void DebyeHuckel::initThermo()
906 {
907  MolalityVPSSTP::initThermo();
908  initLengths();
909 }
910 
911 /*
912  * constructPhaseFile
913  *
914  * Initialization of a Debye-Huckel phase using an
915  * xml file.
916  *
917  * This routine is a precursor to initThermo(XML_Node*)
918  * routine, which does most of the work.
919  *
920  * @param infile XML file containing the description of the
921  * phase
922  *
923  * @param id Optional parameter identifying the name of the
924  * phase. If none is given, the first XML
925  * phase element will be used.
926  */
927 void DebyeHuckel::constructPhaseFile(std::string inputFile, std::string id)
928 {
929 
930  if (inputFile.size() == 0) {
931  throw CanteraError("DebyeHuckel::initThermo",
932  "input file is null");
933  }
934  std::string path = findInputFile(inputFile);
935  ifstream fin(path.c_str());
936  if (!fin) {
937  throw CanteraError("DebyeHuckel::initThermo","could not open "
938  +path+" for reading.");
939  }
940  /*
941  * The phase object automatically constructs an XML object.
942  * Use this object to store information.
943  */
944  XML_Node& phaseNode_XML = xml();
945  XML_Node* fxml = new XML_Node();
946  fxml->build(fin);
947  XML_Node* fxml_phase = findXMLPhase(fxml, id);
948  if (!fxml_phase) {
949  throw CanteraError("DebyeHuckel::initThermo",
950  "ERROR: Can not find phase named " +
951  id + " in file named " + inputFile);
952  }
953  fxml_phase->copy(&phaseNode_XML);
954  constructPhaseXML(*fxml_phase, id);
955  delete fxml;
956 }
957 
958 //! Utility function to assign an integer value from a string for the ElectrolyteSpeciesType field.
959 /*!
960  * @param estString input string that will be interpreted
961  */
962 static int interp_est(std::string estString)
963 {
964  const char* cc = estString.c_str();
965  string lc = lowercase(estString);
966  const char* ccl = lc.c_str();
967  if (!strcmp(ccl, "solvent")) {
968  return cEST_solvent;
969  } else if (!strcmp(ccl, "chargedspecies")) {
970  return cEST_chargedSpecies;
971  } else if (!strcmp(ccl, "weakacidassociated")) {
972  return cEST_weakAcidAssociated;
973  } else if (!strcmp(ccl, "strongacidassociated")) {
974  return cEST_strongAcidAssociated;
975  } else if (!strcmp(ccl, "polarneutral")) {
976  return cEST_polarNeutral;
977  } else if (!strcmp(ccl, "nonpolarneutral")) {
978  return cEST_nonpolarNeutral;
979  }
980  int retn, rval;
981  if ((retn = sscanf(cc, "%d", &rval)) != 1) {
982  return -1;
983  }
984  return rval;
985 }
986 
987 /*
988  * Import and initialize a DebyeHuckel phase
989  * specification in an XML tree into the current object.
990  * Here we read an XML description of the phase.
991  * We import descriptions of the elements that make up the
992  * species in a phase.
993  * We import information about the species, including their
994  * reference state thermodynamic polynomials. We then freeze
995  * the state of the species.
996  *
997  * Then, we read the species molar volumes from the xml
998  * tree to finish the initialization.
999  *
1000  * @param phaseNode This object must be the phase node of a
1001  * complete XML tree
1002  * description of the phase, including all of the
1003  * species data. In other words while "phase" must
1004  * point to an XML phase object, it must have
1005  * sibling nodes "speciesData" that describe
1006  * the species in the phase.
1007  * @param id ID of the phase. If nonnull, a check is done
1008  * to see if phaseNode is pointing to the phase
1009  * with the correct id.
1010  */
1011 void DebyeHuckel::constructPhaseXML(XML_Node& phaseNode, std::string id)
1012 {
1013 
1014  if (id.size() > 0) {
1015  std::string idp = phaseNode.id();
1016  if (idp != id) {
1017  throw CanteraError("DebyeHuckel::constructPhaseXML",
1018  "phasenode and Id are incompatible");
1019  }
1020  }
1021 
1022  /*
1023  * Find the Thermo XML node
1024  */
1025  if (!phaseNode.hasChild("thermo")) {
1026  throw CanteraError("DebyeHuckel::constructPhaseXML",
1027  "no thermo XML node");
1028  }
1029  XML_Node& thermoNode = phaseNode.child("thermo");
1030 
1031  /*
1032  * Possibly change the form of the standard concentrations
1033  */
1034  if (thermoNode.hasChild("standardConc")) {
1035  XML_Node& scNode = thermoNode.child("standardConc");
1036  m_formGC = 2;
1037  std::string formString = scNode.attrib("model");
1038  if (formString != "") {
1039  if (formString == "unity") {
1040  m_formGC = 0;
1041  printf("exit standardConc = unity not done\n");
1042  exit(EXIT_FAILURE);
1043  } else if (formString == "molar_volume") {
1044  m_formGC = 1;
1045  printf("exit standardConc = molar_volume not done\n");
1046  exit(EXIT_FAILURE);
1047  } else if (formString == "solvent_volume") {
1048  m_formGC = 2;
1049  } else {
1050  throw CanteraError("DebyeHuckel::constructPhaseXML",
1051  "Unknown standardConc model: " + formString);
1052  }
1053  }
1054  }
1055  /*
1056  * Get the Name of the Solvent:
1057  * <solvent> solventName </solvent>
1058  */
1059  std::string solventName = "";
1060  if (thermoNode.hasChild("solvent")) {
1061  XML_Node& scNode = thermoNode.child("solvent");
1062  vector<std::string> nameSolventa;
1063  getStringArray(scNode, nameSolventa);
1064  int nsp = static_cast<int>(nameSolventa.size());
1065  if (nsp != 1) {
1066  throw CanteraError("DebyeHuckel::constructPhaseXML",
1067  "badly formed solvent XML node");
1068  }
1069  solventName = nameSolventa[0];
1070  }
1071 
1072  /*
1073  * Determine the form of the Debye-Huckel model,
1074  * m_formDH. We will use this information to size arrays below.
1075  */
1076  if (thermoNode.hasChild("activityCoefficients")) {
1077  XML_Node& scNode = thermoNode.child("activityCoefficients");
1078  m_formDH = DHFORM_DILUTE_LIMIT;
1079  std::string formString = scNode.attrib("model");
1080  if (formString != "") {
1081  if (formString == "Dilute_limit") {
1082  m_formDH = DHFORM_DILUTE_LIMIT;
1083  } else if (formString == "Bdot_with_variable_a") {
1084  m_formDH = DHFORM_BDOT_AK ;
1085  } else if (formString == "Bdot_with_common_a") {
1086  m_formDH = DHFORM_BDOT_ACOMMON;
1087  } else if (formString == "Beta_ij") {
1088  m_formDH = DHFORM_BETAIJ;
1089  } else if (formString == "Pitzer_with_Beta_ij") {
1090  m_formDH = DHFORM_PITZER_BETAIJ;
1091  } else {
1092  throw CanteraError("DebyeHuckel::constructPhaseXML",
1093  "Unknown standardConc model: " + formString);
1094  }
1095  }
1096  } else {
1097  /*
1098  * If there is no XML node named "activityCoefficients", assume
1099  * that we are doing the extreme dilute limit assumption
1100  */
1101  m_formDH = DHFORM_DILUTE_LIMIT;
1102  }
1103 
1104  /*
1105  * Call the Cantera importPhase() function. This will import
1106  * all of the species into the phase. This will also handle
1107  * all of the solvent and solute standard states
1108  */
1109  bool m_ok = importPhase(phaseNode, this);
1110  if (!m_ok) {
1111  throw CanteraError("DebyeHuckel::constructPhaseXML",
1112  "importPhase failed ");
1113  }
1114 }
1115 
1116 /*
1117  * Process the XML file after species are set up.
1118  *
1119  * This gets called from importPhase(). It processes the XML file
1120  * after the species are set up. This is the main routine for
1121  * reading in activity coefficient parameters.
1122  *
1123  * @param phaseNode This object must be the phase node of a
1124  * complete XML tree
1125  * description of the phase, including all of the
1126  * species data. In other words while "phase" must
1127  * point to an XML phase object, it must have
1128  * sibling nodes "speciesData" that describe
1129  * the species in the phase.
1130  * @param id ID of the phase. If nonnull, a check is done
1131  * to see if phaseNode is pointing to the phase
1132  * with the correct id.
1133  */
1134 void DebyeHuckel::
1135 initThermoXML(XML_Node& phaseNode, std::string id)
1136 {
1137  std::string stemp;
1138  /*
1139  * Find the Thermo XML node
1140  */
1141  if (!phaseNode.hasChild("thermo")) {
1142  throw CanteraError("HMWSoln::initThermoXML",
1143  "no thermo XML node");
1144  }
1145  XML_Node& thermoNode = phaseNode.child("thermo");
1146 
1147  /*
1148  * Possibly change the form of the standard concentrations
1149  */
1150  if (thermoNode.hasChild("standardConc")) {
1151  XML_Node& scNode = thermoNode.child("standardConc");
1152  m_formGC = 2;
1153  std::string formString = scNode.attrib("model");
1154  if (formString != "") {
1155  if (formString == "unity") {
1156  m_formGC = 0;
1157  printf("exit standardConc = unity not done\n");
1158  exit(EXIT_FAILURE);
1159  } else if (formString == "molar_volume") {
1160  m_formGC = 1;
1161  printf("exit standardConc = molar_volume not done\n");
1162  exit(EXIT_FAILURE);
1163  } else if (formString == "solvent_volume") {
1164  m_formGC = 2;
1165  } else {
1166  throw CanteraError("DebyeHuckel::constructPhaseXML",
1167  "Unknown standardConc model: " + formString);
1168  }
1169  }
1170  }
1171 
1172 
1173 
1174  /*
1175  * Reconcile the solvent name and index.
1176  */
1177  /*
1178  * Get the Name of the Solvent:
1179  * <solvent> solventName </solvent>
1180  */
1181  std::string solventName = "";
1182  if (thermoNode.hasChild("solvent")) {
1183  XML_Node& scNode = thermoNode.child("solvent");
1184  vector<std::string> nameSolventa;
1185  getStringArray(scNode, nameSolventa);
1186  int nsp = static_cast<int>(nameSolventa.size());
1187  if (nsp != 1) {
1188  throw CanteraError("DebyeHuckel::initThermoXML",
1189  "badly formed solvent XML node");
1190  }
1191  solventName = nameSolventa[0];
1192  }
1193  for (size_t k = 0; k < m_kk; k++) {
1194  std::string sname = speciesName(k);
1195  if (solventName == sname) {
1196  m_indexSolvent = k;
1197  break;
1198  }
1199  }
1200  if (m_indexSolvent == npos) {
1201  cout << "DebyeHuckel::initThermoXML: Solvent Name not found"
1202  << endl;
1203  throw CanteraError("DebyeHuckel::initThermoXML",
1204  "Solvent name not found");
1205  }
1206  if (m_indexSolvent != 0) {
1207  throw CanteraError("DebyeHuckel::initThermoXML",
1208  "Solvent " + solventName +
1209  " should be first species");
1210  }
1211 
1212  /*
1213  * Determine the form of the Debye-Huckel model,
1214  * m_formDH. We will use this information to size arrays below.
1215  */
1216  if (thermoNode.hasChild("activityCoefficients")) {
1217  XML_Node& scNode = thermoNode.child("activityCoefficients");
1218  m_formDH = DHFORM_DILUTE_LIMIT;
1219  std::string formString = scNode.attrib("model");
1220  if (formString != "") {
1221  if (formString == "Dilute_limit") {
1222  m_formDH = DHFORM_DILUTE_LIMIT;
1223  } else if (formString == "Bdot_with_variable_a") {
1224  m_formDH = DHFORM_BDOT_AK ;
1225  } else if (formString == "Bdot_with_common_a") {
1226  m_formDH = DHFORM_BDOT_ACOMMON;
1227  } else if (formString == "Beta_ij") {
1228  m_formDH = DHFORM_BETAIJ;
1229  } else if (formString == "Pitzer_with_Beta_ij") {
1230  m_formDH = DHFORM_PITZER_BETAIJ;
1231  } else {
1232  throw CanteraError("DebyeHuckel::constructPhaseXML",
1233  "Unknown standardConc model: " + formString);
1234  }
1235  }
1236  } else {
1237  /*
1238  * If there is no XML node named "activityCoefficients", assume
1239  * that we are doing the extreme dilute limit assumption
1240  */
1241  m_formDH = DHFORM_DILUTE_LIMIT;
1242  }
1243 
1244  /*
1245  * Initialize all of the lengths of arrays in the object
1246  * now that we know what species are in the phase.
1247  */
1248  initThermo();
1249 
1250  /*
1251  * Now go get the specification of the standard states for
1252  * species in the solution. This includes the molar volumes
1253  * data blocks for incompressible species.
1254  */
1255  XML_Node& speciesList = phaseNode.child("speciesArray");
1256  XML_Node* speciesDB =
1257  get_XML_NameID("speciesData", speciesList["datasrc"],
1258  &phaseNode.root());
1259  const vector<string>&sss = speciesNames();
1260 
1261  for (size_t k = 0; k < m_kk; k++) {
1262  XML_Node* s = speciesDB->findByAttr("name", sss[k]);
1263  if (!s) {
1264  throw CanteraError("DebyeHuckel::initThermoXML",
1265  "Species Data Base " + sss[k] + " not found");
1266  }
1267  XML_Node* ss = s->findByName("standardState");
1268  if (!ss) {
1269  throw CanteraError("DebyeHuckel::initThermoXML",
1270  "Species " + sss[k] +
1271  " standardState XML block not found");
1272  }
1273  std::string modelStringa = ss->attrib("model");
1274  if (modelStringa == "") {
1275  throw CanteraError("DebyeHuckel::initThermoXML",
1276  "Species " + sss[k] +
1277  " standardState XML block model attribute not found");
1278  }
1279  std::string modelString = lowercase(modelStringa);
1280 
1281  if (k == 0) {
1282  if (modelString == "wateriapws" || modelString == "real_water" ||
1283  modelString == "waterpdss") {
1284  /*
1285  * Initialize the water standard state model
1286  */
1287  m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0)) ;
1288  if (!m_waterSS) {
1289  throw CanteraError("HMWSoln::installThermoXML",
1290  "Dynamic cast to PDSS_Water failed");
1291  }
1292  /*
1293  * Fill in the molar volume of water (m3/kmol)
1294  * at standard conditions to fill in the m_speciesSize entry
1295  * with something reasonable.
1296  */
1297  m_waterSS->setState_TP(300., OneAtm);
1298  double dens = m_waterSS->density();
1299  double mw = m_waterSS->molecularWeight();
1300  m_speciesSize[0] = mw / dens;
1301 #ifdef DEBUG_MODE_NOT
1302  cout << "Solvent species " << sss[k] << " has volume " <<
1303  m_speciesSize[k] << endl;
1304 #endif
1305  } else if (modelString == "constant_incompressible") {
1306  m_speciesSize[k] = getFloat(*ss, "molarVolume", "toSi");
1307 #ifdef DEBUG_MODE_NOT
1308  cout << "species " << sss[k] << " has volume " <<
1309  m_speciesSize[k] << endl;
1310 #endif
1311  } else {
1312  throw CanteraError("DebyeHuckel::initThermoXML",
1313  "Solvent SS Model \"" + modelStringa +
1314  "\" is not known");
1315  }
1316  } else {
1317  if (modelString != "constant_incompressible") {
1318  throw CanteraError("DebyeHuckel::initThermoXML",
1319  "Solute SS Model \"" + modelStringa +
1320  "\" is not known");
1321  }
1322  m_speciesSize[k] = getFloat(*ss, "molarVolume", "toSI");
1323 #ifdef DEBUG_MODE_NOT
1324  cout << "species " << sss[k] << " has volume " <<
1325  m_speciesSize[k] << endl;
1326 #endif
1327  }
1328 
1329  }
1330 
1331  /*
1332  * Go get all of the coefficients and factors in the
1333  * activityCoefficients XML block
1334  */
1335  XML_Node* acNodePtr = 0;
1336  if (thermoNode.hasChild("activityCoefficients")) {
1337  XML_Node& acNode = thermoNode.child("activityCoefficients");
1338  acNodePtr = &acNode;
1339  /*
1340  * Look for parameters for A_Debye
1341  */
1342  if (acNode.hasChild("A_Debye")) {
1343  XML_Node* ss = acNode.findByName("A_Debye");
1344  string modelStringa = ss->attrib("model");
1345  string modelString = lowercase(modelStringa);
1346  if (modelString != "") {
1347  if (modelString == "water") {
1348  m_form_A_Debye = A_DEBYE_WATER;
1349  } else {
1350  throw CanteraError("DebyeHuckel::initThermoXML",
1351  "A_Debye Model \"" + modelStringa +
1352  "\" is not known");
1353  }
1354  } else {
1355  m_A_Debye = getFloat(acNode, "A_Debye");
1356 #ifdef DEBUG_HKM_NOT
1357  cout << "A_Debye = " << m_A_Debye << endl;
1358 #endif
1359  }
1360  }
1361 
1362  /*
1363  * Initialize the water property calculator. It will share
1364  * the internal eos water calculator.
1365  */
1366  if (m_form_A_Debye == A_DEBYE_WATER) {
1367  if (m_waterProps) {
1368  delete m_waterProps;
1369  }
1370  m_waterProps = new WaterProps(m_waterSS);
1371  }
1372 
1373  /*
1374  * Look for parameters for B_Debye
1375  */
1376  if (acNode.hasChild("B_Debye")) {
1377  m_B_Debye = getFloat(acNode, "B_Debye");
1378 #ifdef DEBUG_HKM_NOT
1379  cout << "B_Debye = " << m_B_Debye << endl;
1380 #endif
1381  }
1382 
1383  /*
1384  * Look for parameters for B_dot
1385  */
1386  if (acNode.hasChild("B_dot")) {
1387  if (m_formDH == DHFORM_BETAIJ ||
1388  m_formDH == DHFORM_DILUTE_LIMIT ||
1389  m_formDH == DHFORM_PITZER_BETAIJ) {
1390  throw CanteraError("DebyeHuckel:init",
1391  "B_dot entry in the wrong DH form");
1392  }
1393  double bdot_common = getFloat(acNode, "B_dot");
1394 #ifdef DEBUG_HKM_NOT
1395  cout << "B_dot = " << bdot_common << endl;
1396 #endif
1397  /*
1398  * Set B_dot parameters for charged species
1399  */
1400  for (size_t k = 0; k < m_kk; k++) {
1401  double z_k = charge(k);
1402  if (fabs(z_k) > 0.0001) {
1403  m_B_Dot[k] = bdot_common;
1404  } else {
1405  m_B_Dot[k] = 0.0;
1406  }
1407  }
1408  }
1409 
1410  /*
1411  * Look for Parameters for the Maximum Ionic Strength
1412  */
1413  if (acNode.hasChild("maxIonicStrength")) {
1414  m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
1415 #ifdef DEBUG_HKM_NOT
1416  cout << "m_maxIionicStrength = "
1417  <<m_maxIionicStrength << endl;
1418 #endif
1419  }
1420 
1421  /*
1422  * Look for Helgeson Parameters
1423  */
1424  if (acNode.hasChild("UseHelgesonFixedForm")) {
1425  m_useHelgesonFixedForm = true;
1426  } else {
1427  m_useHelgesonFixedForm = false;
1428  }
1429 
1430  /*
1431  * Look for parameters for the Ionic radius
1432  */
1433  if (acNode.hasChild("ionicRadius")) {
1434  XML_Node& irNode = acNode.child("ionicRadius");
1435 
1436  std::string Aunits = "";
1437  double Afactor = 1.0;
1438  if (irNode.hasAttrib("units")) {
1439  std::string Aunits = irNode.attrib("units");
1440  Afactor = toSI(Aunits);
1441  }
1442 
1443  if (irNode.hasAttrib("default")) {
1444  std::string ads = irNode.attrib("default");
1445  double ad = fpValue(ads);
1446  for (size_t k = 0; k < m_kk; k++) {
1447  m_Aionic[k] = ad * Afactor;
1448  }
1449  }
1450 
1451  /*
1452  * If the Debye-Huckel form is BDOT_AK, we can
1453  * have separate values for the denominator's ionic
1454  * size. -> That's how the activity coefficient is
1455  * parameterized. In this case only do we allow the
1456  * code to read in these parameters.
1457  */
1458  if (m_formDH == DHFORM_BDOT_AK) {
1459  /*
1460  * Define a string-string map, and interpret the
1461  * value of the xml element as binary pairs separated
1462  * by colons, e.g.:
1463  * Na+:3.0
1464  * Cl-:4.0
1465  * H+:9.0
1466  * OH-:3.5
1467  * Read them into the map.
1468  */
1469  map<string, string> m;
1470  getMap(irNode, m);
1471  /*
1472  * Iterate over the map pairs, interpreting the
1473  * first string as a species in the current phase.
1474  * If no match is made, silently ignore the
1475  * lack of agreement (HKM -> may be changed in the
1476  * future).
1477  */
1478  map<std::string,std::string>::const_iterator _b = m.begin();
1479  for (; _b != m.end(); ++_b) {
1480  size_t kk = speciesIndex(_b->first);
1481  m_Aionic[kk] = fpValue(_b->second) * Afactor;
1482 
1483  }
1484  }
1485  }
1486  /*
1487  * Get the matrix of coefficients for the Beta
1488  * binary interaction parameters. We assume here that
1489  * this matrix is symmetric, so that we only have to
1490  * input 1/2 of the values.
1491  */
1492  if (acNode.hasChild("DHBetaMatrix")) {
1493  if (m_formDH == DHFORM_BETAIJ ||
1494  m_formDH == DHFORM_PITZER_BETAIJ) {
1495  XML_Node& irNode = acNode.child("DHBetaMatrix");
1496  const vector<string>& sn = speciesNames();
1497  getMatrixValues(irNode, sn, sn, m_Beta_ij, true, true);
1498  } else {
1499  throw CanteraError("DebyeHuckel::initThermoXML:",
1500  "DHBetaMatrix found for wrong type");
1501  }
1502  }
1503 
1504  /*
1505  * Fill in parameters for the calculation of the
1506  * stoichiometric Ionic Strength
1507  *
1508  * The default is that stoich charge is the same as the
1509  * regular charge.
1510  */
1511  m_speciesCharge_Stoich.resize(m_kk, 0.0);
1512  for (size_t k = 0; k < m_kk; k++) {
1513  m_speciesCharge_Stoich[k] = m_speciesCharge[k];
1514  }
1515  /*
1516  * First look at the species database.
1517  * -> Look for the subelement "stoichIsMods"
1518  * in each of the species SS databases.
1519  */
1520  std::vector<const XML_Node*> xspecies= speciesData();
1521  std::string kname, jname;
1522  size_t jj = xspecies.size();
1523  for (size_t k = 0; k < m_kk; k++) {
1524  size_t jmap = npos;
1525  kname = speciesName(k);
1526  for (size_t j = 0; j < jj; j++) {
1527  const XML_Node& sp = *xspecies[j];
1528  jname = sp["name"];
1529  if (jname == kname) {
1530  jmap = j;
1531  break;
1532  }
1533  }
1534  if (jmap != npos) {
1535  const XML_Node& sp = *xspecies[jmap];
1536  if (sp.hasChild("stoichIsMods")) {
1537  double val = getFloat(sp, "stoichIsMods");
1538  m_speciesCharge_Stoich[k] = val;
1539  }
1540  }
1541  }
1542 
1543  /*
1544  * Now look at the activity coefficient database
1545  */
1546  if (acNodePtr) {
1547  if (acNodePtr->hasChild("stoichIsMods")) {
1548  XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
1549 
1550  map<std::string, std::string> msIs;
1551  getMap(sIsNode, msIs);
1552  map<std::string,std::string>::const_iterator _b = msIs.begin();
1553  for (; _b != msIs.end(); ++_b) {
1554  size_t kk = speciesIndex(_b->first);
1555  double val = fpValue(_b->second);
1556  m_speciesCharge_Stoich[kk] = val;
1557  }
1558  }
1559  }
1560  }
1561 
1562  /*
1563  * Fill in the vector specifying the electrolyte species
1564  * type
1565  *
1566  * First fill in default values. Everthing is either
1567  * a charge species, a nonpolar neutral, or the solvent.
1568  */
1569  for (size_t k = 0; k < m_kk; k++) {
1570  if (fabs(m_speciesCharge[k]) > 0.0001) {
1571  m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
1572  if (fabs(m_speciesCharge_Stoich[k] - m_speciesCharge[k])
1573  > 0.0001) {
1574  m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
1575  }
1576  } else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) {
1577  m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
1578  } else {
1579  m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
1580  }
1581  }
1582  m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
1583  /*
1584  * First look at the species database.
1585  * -> Look for the subelement "stoichIsMods"
1586  * in each of the species SS databases.
1587  */
1588  std::vector<const XML_Node*> xspecies= speciesData();
1589  const XML_Node* spPtr = 0;
1590  std::string kname;
1591  for (size_t k = 0; k < m_kk; k++) {
1592  kname = speciesName(k);
1593  spPtr = xspecies[k];
1594  if (!spPtr) {
1595  if (spPtr->hasChild("electrolyteSpeciesType")) {
1596  std::string est = getChildValue(*spPtr, "electrolyteSpeciesType");
1597  if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
1598  throw CanteraError("DebyeHuckel:initThermoXML",
1599  "Bad electrolyte type: " + est);
1600  }
1601  }
1602  }
1603  }
1604  /*
1605  * Then look at the phase thermo specification
1606  */
1607  if (acNodePtr) {
1608  if (acNodePtr->hasChild("electrolyteSpeciesType")) {
1609  XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
1610  map<std::string, std::string> msEST;
1611  getMap(ESTNode, msEST);
1612  map<std::string,std::string>::const_iterator _b = msEST.begin();
1613  for (; _b != msEST.end(); ++_b) {
1614  size_t kk = speciesIndex(_b->first);
1615  std::string est = _b->second;
1616  if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
1617  throw CanteraError("DebyeHuckel:initThermoXML",
1618  "Bad electrolyte type: " + est);
1619  }
1620  }
1621  }
1622  }
1623 
1624 
1625 
1626  /*
1627  * Lastly set the state
1628  */
1629  if (phaseNode.hasChild("state")) {
1630  XML_Node& stateNode = phaseNode.child("state");
1631  setStateFromXML(stateNode);
1632  }
1633 
1634 }
1635 
1636 /*
1637  * @internal
1638  * Set equation of state parameters. The number and meaning of
1639  * these depends on the subclass.
1640  * @param n number of parameters
1641  * @param c array of \i n coefficients
1642  *
1643  */
1644 void DebyeHuckel::setParameters(int n, doublereal* const c)
1645 {
1646 }
1647 
1648 void DebyeHuckel::getParameters(int& n, doublereal* const c) const
1649 {
1650 }
1651 
1652 /*
1653  * Set equation of state parameter values from XML
1654  * entries. This method is called by function importPhase in
1655  * file importCTML.cpp when processing a phase definition in
1656  * an input file. It should be overloaded in subclasses to set
1657  * any parameters that are specific to that particular phase
1658  * model.
1659  *
1660  * @param eosdata An XML_Node object corresponding to
1661  * the "thermo" entry for this phase in the input file.
1662  *
1663  * HKM -> Right now, the parameters are set elsewhere (initThermoXML)
1664  * It just didn't seem to fit.
1665  */
1666 void DebyeHuckel::setParametersFromXML(const XML_Node& eosdata)
1667 {
1668 }
1669 
1670 /*
1671  * Report the molar volume of species k
1672  *
1673  * units - \f$ m^3 kmol^-1 \f$
1674  */
1675 // double DebyeHuckel::speciesMolarVolume(int k) const {
1676 // return m_speciesSize[k];
1677 //}
1678 
1679 
1680 /*
1681  * A_Debye_TP() (virtual)
1682  *
1683  * Returns the A_Debye parameter as a function of temperature
1684  * and pressure.
1685  *
1686  * The default is to assume that it is constant, given
1687  * in the initialization process and stored in the
1688  * member double, m_A_Debye
1689  */
1690 double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const
1691 {
1692  double T = temperature();
1693  double A;
1694  if (tempArg != -1.0) {
1695  T = tempArg;
1696  }
1697  double P = pressure();
1698  if (presArg != -1.0) {
1699  P = presArg;
1700  }
1701 
1702  switch (m_form_A_Debye) {
1703  case A_DEBYE_CONST:
1704  A = m_A_Debye;
1705  break;
1706  case A_DEBYE_WATER:
1707  A = m_waterProps->ADebye(T, P, 0);
1708  m_A_Debye = A;
1709  break;
1710  default:
1711  printf("shouldn't be here\n");
1712  exit(EXIT_FAILURE);
1713  }
1714  return A;
1715 }
1716 
1717 /*
1718  * dA_DebyedT_TP() (virtual)
1719  *
1720  * Returns the derivative of the A_Debye parameter with
1721  * respect to temperature as a function of temperature
1722  * and pressure.
1723  *
1724  * units = A_Debye has units of sqrt(gmol kg-1).
1725  * Temp has units of Kelvin.
1726  */
1727 double DebyeHuckel::dA_DebyedT_TP(double tempArg, double presArg) const
1728 {
1729  double T = temperature();
1730  if (tempArg != -1.0) {
1731  T = tempArg;
1732  }
1733  double P = pressure();
1734  if (presArg != -1.0) {
1735  P = presArg;
1736  }
1737  double dAdT;
1738  switch (m_form_A_Debye) {
1739  case A_DEBYE_CONST:
1740  dAdT = 0.0;
1741  break;
1742  case A_DEBYE_WATER:
1743  dAdT = m_waterProps->ADebye(T, P, 1);
1744  break;
1745  default:
1746  printf("shouldn't be here\n");
1747  exit(EXIT_FAILURE);
1748  }
1749  return dAdT;
1750 }
1751 
1752 /*
1753  * d2A_DebyedT2_TP() (virtual)
1754  *
1755  * Returns the 2nd derivative of the A_Debye parameter with
1756  * respect to temperature as a function of temperature
1757  * and pressure.
1758  *
1759  * units = A_Debye has units of sqrt(gmol kg-1).
1760  * Temp has units of Kelvin.
1761  */
1762 double DebyeHuckel::d2A_DebyedT2_TP(double tempArg, double presArg) const
1763 {
1764  double T = temperature();
1765  if (tempArg != -1.0) {
1766  T = tempArg;
1767  }
1768  double P = pressure();
1769  if (presArg != -1.0) {
1770  P = presArg;
1771  }
1772  double d2AdT2;
1773  switch (m_form_A_Debye) {
1774  case A_DEBYE_CONST:
1775  d2AdT2 = 0.0;
1776  break;
1777  case A_DEBYE_WATER:
1778  d2AdT2 = m_waterProps->ADebye(T, P, 2);
1779  break;
1780  default:
1781  printf("shouldn't be here\n");
1782  exit(EXIT_FAILURE);
1783  }
1784  return d2AdT2;
1785 }
1786 
1787 /*
1788  * dA_DebyedP_TP() (virtual)
1789  *
1790  * Returns the derivative of the A_Debye parameter with
1791  * respect to pressure, as a function of temperature
1792  * and pressure.
1793  *
1794  * units = A_Debye has units of sqrt(gmol kg-1).
1795  * Pressure has units of pascals.
1796  */
1797 double DebyeHuckel::dA_DebyedP_TP(double tempArg, double presArg) const
1798 {
1799  double T = temperature();
1800  if (tempArg != -1.0) {
1801  T = tempArg;
1802  }
1803  double P = pressure();
1804  if (presArg != -1.0) {
1805  P = presArg;
1806  }
1807  double dAdP;
1808  switch (m_form_A_Debye) {
1809  case A_DEBYE_CONST:
1810  dAdP = 0.0;
1811  break;
1812  case A_DEBYE_WATER:
1813  dAdP = m_waterProps->ADebye(T, P, 3);
1814  break;
1815  default:
1816  printf("shouldn't be here\n");
1817  exit(EXIT_FAILURE);
1818  }
1819  return dAdP;
1820 }
1821 
1822 /*
1823  * ----------- Critical State Properties --------------------------
1824  */
1825 
1826 /*
1827  * ---------- Other Property Functions
1828  */
1829 double DebyeHuckel::AionicRadius(int k) const
1830 {
1831  return m_Aionic[k];
1832 }
1833 
1834 /*
1835  * ------------ Private and Restricted Functions ------------------
1836  */
1837 
1838 /*
1839  * Bail out of functions with an error exit if they are not
1840  * implemented.
1841  */
1842 doublereal DebyeHuckel::err(std::string msg) const
1843 {
1844  throw CanteraError("DebyeHuckel",
1845  "Unfinished func called: " + msg);
1846  return 0.0;
1847 }
1848 
1849 
1850 /*
1851  * initLengths():
1852  *
1853  * This internal function adjusts the lengths of arrays based on
1854  * the number of species
1855  */
1856 void DebyeHuckel::initLengths()
1857 {
1858  m_kk = nSpecies();
1859 
1860  /*
1861  * Obtain the limits of the temperature from the species
1862  * thermo handler's limits.
1863  */
1864  m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
1865  m_speciesSize.resize(m_kk);
1866  m_Aionic.resize(m_kk, 0.0);
1867  m_lnActCoeffMolal.resize(m_kk, 0.0);
1868  m_dlnActCoeffMolaldT.resize(m_kk, 0.0);
1869  m_d2lnActCoeffMolaldT2.resize(m_kk, 0.0);
1870  m_dlnActCoeffMolaldP.resize(m_kk, 0.0);
1871  m_B_Dot.resize(m_kk, 0.0);
1872  m_expg0_RT.resize(m_kk, 0.0);
1873  m_pe.resize(m_kk, 0.0);
1874  m_pp.resize(m_kk, 0.0);
1875  m_tmpV.resize(m_kk, 0.0);
1876  if (m_formDH == DHFORM_BETAIJ ||
1877  m_formDH == DHFORM_PITZER_BETAIJ) {
1878  m_Beta_ij.resize(m_kk, m_kk, 0.0);
1879  }
1880 }
1881 
1882 /*
1883  * nonpolarActCoeff() (private)
1884  *
1885  * Static function that implements the non-polar species
1886  * salt-out modifications.
1887  * Returns the calculated activity coefficients.
1888  */
1889 double DebyeHuckel::_nonpolarActCoeff(double IionicMolality) const
1890 {
1891  double I2 = IionicMolality * IionicMolality;
1892  double l10actCoeff =
1893  m_npActCoeff[0] * IionicMolality +
1894  m_npActCoeff[1] * I2 +
1895  m_npActCoeff[2] * I2 * IionicMolality;
1896  return pow(10.0 , l10actCoeff);
1897 }
1898 
1899 
1900 /**
1901  * _osmoticCoeffHelgesonFixedForm()
1902  *
1903  * Formula for the osmotic coefficient that occurs in
1904  * the GWB. It is originally from Helgeson for a variable
1905  * NaCl brine. It's to be used with extreme caution.
1906  */
1907 double DebyeHuckel::
1908 _osmoticCoeffHelgesonFixedForm() const
1909 {
1910  const double a0 = 1.454;
1911  const double b0 = 0.02236;
1912  const double c0 = 9.380E-3;
1913  const double d0 = -5.362E-4;
1914  double Is = m_IionicMolalityStoich;
1915  if (Is <= 0.0) {
1916  return 0.0;
1917  }
1918  double Is2 = Is * Is;
1919  double bhat = 1.0 + a0 * sqrt(Is);
1920  double func = bhat - 2.0 * log(bhat) - 1.0/bhat;
1921  double v1 = m_A_Debye / (a0 * a0 * a0 * Is) * func;
1922  double oc = 1.0 - v1 + b0 * Is / 2.0 + 2.0 * c0 * Is2 / 3.0
1923  + 3.0 * d0 * Is2 * Is / 4.0;
1924  return oc;
1925 }
1926 
1927 
1928 /*
1929  * _activityWaterHelgesonFixedForm()
1930  *
1931  * Formula for the log of the activity of the water
1932  * solvent that occurs in
1933  * the GWB. It is originally from Helgeson for a variable
1934  * NaCl brine. It's to be used with extreme caution.
1935  */
1936 double DebyeHuckel::
1937 _lnactivityWaterHelgesonFixedForm() const
1938 {
1939  /*
1940  * Update the internally stored vector of molalities
1941  */
1942  calcMolalities();
1943  double oc = _osmoticCoeffHelgesonFixedForm();
1944  double sum = 0.0;
1945  for (size_t k = 0; k < m_kk; k++) {
1946  if (k != m_indexSolvent) {
1947  sum += std::max(m_molalities[k], 0.0);
1948  }
1949  }
1950  if (sum > 2.0 * m_maxIionicStrength) {
1951  sum = 2.0 * m_maxIionicStrength;
1952  };
1953  double lac = - m_Mnaught * sum * oc;
1954  return lac;
1955 }
1956 
1957 /*
1958  * s_update_lnMolalityActCoeff():
1959  *
1960  * Using internally stored values, this function calculates
1961  * the activity coefficients for all species.
1962  *
1963  * The ln(activity_solvent) is first calculated for the
1964  * solvent. Then the molar based activity coefficient
1965  * is calculated and returned.
1966  *
1967  * ( Note this is the main routine for implementing the
1968  * activity coefficient formulation.)
1969  */
1970 void DebyeHuckel::s_update_lnMolalityActCoeff() const
1971 {
1972  double z_k, zs_k1, zs_k2;
1973  /*
1974  * Update the internally stored vector of molalities
1975  */
1976  calcMolalities();
1977  /*
1978  * Calculate the apparent (real) ionic strength.
1979  *
1980  * Note this is not the stoichiometric ionic strengh,
1981  * where reactions of ions forming neutral salts
1982  * are ignorred in calculating the ionic strength.
1983  */
1984  m_IionicMolality = 0.0;
1985  for (size_t k = 0; k < m_kk; k++) {
1986  z_k = m_speciesCharge[k];
1987  m_IionicMolality += m_molalities[k] * z_k * z_k;
1988  }
1989  m_IionicMolality /= 2.0;
1990 
1991  if (m_IionicMolality > m_maxIionicStrength) {
1992  m_IionicMolality = m_maxIionicStrength;
1993  }
1994 
1995  /*
1996  * Calculate the stoichiometric ionic charge
1997  */
1998  m_IionicMolalityStoich = 0.0;
1999  for (size_t k = 0; k < m_kk; k++) {
2000  z_k = m_speciesCharge[k];
2001  zs_k1 = m_speciesCharge_Stoich[k];
2002  if (z_k == zs_k1) {
2003  m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
2004  } else {
2005  zs_k2 = z_k - zs_k1;
2006  m_IionicMolalityStoich
2007  += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
2008  }
2009  }
2010  m_IionicMolalityStoich /= 2.0;
2011 
2012  if (m_IionicMolalityStoich > m_maxIionicStrength) {
2013  m_IionicMolalityStoich = m_maxIionicStrength;
2014  }
2015 
2016  /*
2017  * Possibly update the stored value of the
2018  * Debye-Huckel parameter A_Debye
2019  * This parameter appears on the top of the activity
2020  * coefficient expression.
2021  * It depends on T (and P), as it depends explicity
2022  * on the temperature. Also, the dielectric constant
2023  * is usually a fairly strong function of T, also.
2024  */
2025  m_A_Debye = A_Debye_TP();
2026 
2027  /*
2028  * Calculate a safe value for the mole fraction
2029  * of the solvent
2030  */
2031  double xmolSolvent = moleFraction(m_indexSolvent);
2032  xmolSolvent = std::max(8.689E-3, xmolSolvent);
2033 
2034  int est;
2035  double ac_nonPolar = 1.0;
2036  double numTmp = m_A_Debye * sqrt(m_IionicMolality);
2037  double denomTmp = m_B_Debye * sqrt(m_IionicMolality);
2038  double coeff;
2039  double lnActivitySolvent = 0.0;
2040  double tmp;
2041  double tmpLn;
2042  double y, yp1, sigma;
2043  switch (m_formDH) {
2044  case DHFORM_DILUTE_LIMIT:
2045  for (size_t k = 0; k < m_kk; k++) {
2046  z_k = m_speciesCharge[k];
2047  m_lnActCoeffMolal[k] = - z_k * z_k * numTmp;
2048  }
2049  lnActivitySolvent =
2050  (xmolSolvent - 1.0)/xmolSolvent +
2051  2.0 / 3.0 * m_A_Debye * m_Mnaught *
2052  m_IionicMolality * sqrt(m_IionicMolality);
2053  break;
2054 
2055  case DHFORM_BDOT_AK:
2056  ac_nonPolar = _nonpolarActCoeff(m_IionicMolality);
2057  for (size_t k = 0; k < m_kk; k++) {
2058  est = m_electrolyteSpeciesType[k];
2059  if (est == cEST_nonpolarNeutral) {
2060  m_lnActCoeffMolal[k] = log(ac_nonPolar);
2061  } else {
2062  z_k = m_speciesCharge[k];
2063  m_lnActCoeffMolal[k] =
2064  - z_k * z_k * numTmp / (1.0 + denomTmp * m_Aionic[k])
2065  + log(10.0) * m_B_Dot[k] * m_IionicMolality;
2066  }
2067  }
2068 
2069  lnActivitySolvent = (xmolSolvent - 1.0)/xmolSolvent;
2070  coeff = 2.0 / 3.0 * m_A_Debye * m_Mnaught
2071  * sqrt(m_IionicMolality);
2072  tmp = 0.0;
2073  if (denomTmp > 0.0) {
2074  for (size_t k = 0; k < m_kk; k++) {
2075  if (k != m_indexSolvent || m_Aionic[k] != 0.0) {
2076  y = denomTmp * m_Aionic[k];
2077  yp1 = y + 1.0;
2078  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2079  z_k = m_speciesCharge[k];
2080  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
2081  }
2082  }
2083  }
2084  lnActivitySolvent += coeff * tmp;
2085  tmp = 0.0;
2086  for (size_t k = 0; k < m_kk; k++) {
2087  z_k = m_speciesCharge[k];
2088  if ((k != m_indexSolvent) && (z_k != 0.0)) {
2089  tmp += m_B_Dot[k] * m_molalities[k];
2090  }
2091  }
2092  lnActivitySolvent -=
2093  m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
2094 
2095  /*
2096  * Special section to implement the Helgeson fixed form
2097  * for the water brine activity coefficient.
2098  */
2099  if (m_useHelgesonFixedForm) {
2100  lnActivitySolvent = _lnactivityWaterHelgesonFixedForm();
2101  }
2102  break;
2103 
2104  case DHFORM_BDOT_ACOMMON:
2105  denomTmp *= m_Aionic[0];
2106  for (size_t k = 0; k < m_kk; k++) {
2107  z_k = m_speciesCharge[k];
2108  m_lnActCoeffMolal[k] =
2109  - z_k * z_k * numTmp / (1.0 + denomTmp)
2110  + log(10.0) * m_B_Dot[k] * m_IionicMolality;
2111  }
2112  if (denomTmp > 0.0) {
2113  y = denomTmp;
2114  yp1 = y + 1.0;
2115  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2116  } else {
2117  sigma = 0.0;
2118  }
2119  lnActivitySolvent =
2120  (xmolSolvent - 1.0)/xmolSolvent +
2121  2.0 /3.0 * m_A_Debye * m_Mnaught *
2122  m_IionicMolality * sqrt(m_IionicMolality) * sigma;
2123  tmp = 0.0;
2124  for (size_t k = 0; k < m_kk; k++) {
2125  z_k = m_speciesCharge[k];
2126  if ((k != m_indexSolvent) && (z_k != 0.0)) {
2127  tmp += m_B_Dot[k] * m_molalities[k];
2128  }
2129  }
2130  lnActivitySolvent -=
2131  m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
2132 
2133  break;
2134 
2135  case DHFORM_BETAIJ:
2136  denomTmp = m_B_Debye * m_Aionic[0];
2137  denomTmp *= sqrt(m_IionicMolality);
2138  lnActivitySolvent =
2139  (xmolSolvent - 1.0)/xmolSolvent;
2140 
2141  for (size_t k = 0; k < m_kk; k++) {
2142  if (k != m_indexSolvent) {
2143  z_k = m_speciesCharge[k];
2144  m_lnActCoeffMolal[k] =
2145  - z_k * z_k * numTmp / (1.0 + denomTmp);
2146  for (size_t j = 0; j < m_kk; j++) {
2147  double beta = m_Beta_ij.value(k, j);
2148 #ifdef DEBUG_HKM_NOT
2149  if (beta != 0.0) {
2150  printf("b: k = %d, j = %d, betakj = %g\n",
2151  k, j, beta);
2152  }
2153 #endif
2154  m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
2155  }
2156  }
2157  }
2158  if (denomTmp > 0.0) {
2159  y = denomTmp;
2160  yp1 = y + 1.0;
2161  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1));
2162  } else {
2163  sigma = 0.0;
2164  }
2165  lnActivitySolvent =
2166  (xmolSolvent - 1.0)/xmolSolvent +
2167  2.0 /3.0 * m_A_Debye * m_Mnaught *
2168  m_IionicMolality * sqrt(m_IionicMolality) * sigma;
2169  tmp = 0.0;
2170  for (size_t k = 0; k < m_kk; k++) {
2171  for (size_t j = 0; j < m_kk; j++) {
2172  tmp +=
2173  m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
2174  }
2175  }
2176  lnActivitySolvent -= m_Mnaught * tmp;
2177  break;
2178 
2179  case DHFORM_PITZER_BETAIJ:
2180  denomTmp = m_B_Debye * sqrt(m_IionicMolality);
2181  denomTmp *= m_Aionic[0];
2182  numTmp = m_A_Debye * sqrt(m_IionicMolality);
2183  tmpLn = log(1.0 + denomTmp);
2184  for (size_t k = 0; k < m_kk; k++) {
2185  if (k != m_indexSolvent) {
2186  z_k = m_speciesCharge[k];
2187  m_lnActCoeffMolal[k] =
2188  - z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
2189  m_lnActCoeffMolal[k] +=
2190  - 2.0 * z_k * z_k * m_A_Debye * tmpLn /
2191  (3.0 * m_B_Debye * m_Aionic[0]);
2192  for (size_t j = 0; j < m_kk; j++) {
2193  m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
2194  m_Beta_ij.value(k, j);
2195  }
2196  }
2197  }
2198  sigma = 1.0 / (1.0 + denomTmp);
2199  lnActivitySolvent =
2200  (xmolSolvent - 1.0)/xmolSolvent +
2201  2.0 /3.0 * m_A_Debye * m_Mnaught *
2202  m_IionicMolality * sqrt(m_IionicMolality) * sigma;
2203  tmp = 0.0;
2204  for (size_t k = 0; k < m_kk; k++) {
2205  for (size_t j = 0; j < m_kk; j++) {
2206  tmp +=
2207  m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
2208  }
2209  }
2210  lnActivitySolvent -= m_Mnaught * tmp;
2211  break;
2212 
2213  default:
2214  printf("ERROR\n");
2215  exit(EXIT_FAILURE);
2216  }
2217  /*
2218  * Above, we calculated the ln(activitySolvent). Translate that
2219  * into the molar-based activity coefficient by dividing by
2220  * the solvent mole fraction. Solvents are not on the molality
2221  * scale.
2222  */
2223  xmolSolvent = moleFraction(m_indexSolvent);
2224  m_lnActCoeffMolal[m_indexSolvent] =
2225  lnActivitySolvent - log(xmolSolvent);
2226 }
2227 
2228 /*
2229  * s_update_dMolalityActCoeff_dT() (private, const )
2230  *
2231  * Using internally stored values, this function calculates
2232  * the temperature derivative of the logarithm of the
2233  * activity coefficient for all species in the mechanism.
2234  *
2235  * We assume that the activity coefficients are current.
2236  *
2237  * solvent activity coefficient is on the molality
2238  * scale. Its derivative is too.
2239  */
2240 void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
2241 {
2242  double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
2243  // First we store dAdT explicitly here
2244  double dAdT = dA_DebyedT_TP();
2245  if (dAdT == 0.0) {
2246  for (size_t k = 0; k < m_kk; k++) {
2247  m_dlnActCoeffMolaldT[k] = 0.0;
2248  }
2249  return;
2250  }
2251  /*
2252  * Calculate a safe value for the mole fraction
2253  * of the solvent
2254  */
2255  double xmolSolvent = moleFraction(m_indexSolvent);
2256  xmolSolvent = std::max(8.689E-3, xmolSolvent);
2257 
2258 
2259  double sqrtI = sqrt(m_IionicMolality);
2260  double numdAdTTmp = dAdT * sqrtI;
2261  double denomTmp = m_B_Debye * sqrtI;
2262  double d_lnActivitySolvent_dT = 0;
2263 
2264  switch (m_formDH) {
2265  case DHFORM_DILUTE_LIMIT:
2266  for (size_t k = 1; k < m_kk; k++) {
2267  m_dlnActCoeffMolaldT[k] =
2268  m_lnActCoeffMolal[k] * dAdT / m_A_Debye;
2269  }
2270  d_lnActivitySolvent_dT = 2.0 / 3.0 * dAdT * m_Mnaught *
2271  m_IionicMolality * sqrt(m_IionicMolality);
2272  m_dlnActCoeffMolaldT[m_indexSolvent] = d_lnActivitySolvent_dT;
2273  break;
2274 
2275  case DHFORM_BDOT_AK:
2276  for (size_t k = 0; k < m_kk; k++) {
2277  z_k = m_speciesCharge[k];
2278  m_dlnActCoeffMolaldT[k] =
2279  - z_k * z_k * numdAdTTmp / (1.0 + denomTmp * m_Aionic[k]);
2280  }
2281 
2282  m_dlnActCoeffMolaldT[m_indexSolvent] = 0.0;
2283 
2284  coeff = 2.0 / 3.0 * dAdT * m_Mnaught * sqrtI;
2285  tmp = 0.0;
2286  if (denomTmp > 0.0) {
2287  for (size_t k = 0; k < m_kk; k++) {
2288  y = denomTmp * m_Aionic[k];
2289  yp1 = y + 1.0;
2290  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2291  z_k = m_speciesCharge[k];
2292  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
2293  }
2294  }
2295  m_dlnActCoeffMolaldT[m_indexSolvent] += coeff * tmp;
2296  break;
2297 
2298  case DHFORM_BDOT_ACOMMON:
2299  denomTmp *= m_Aionic[0];
2300  for (size_t k = 0; k < m_kk; k++) {
2301  z_k = m_speciesCharge[k];
2302  m_dlnActCoeffMolaldT[k] =
2303  - z_k * z_k * numdAdTTmp / (1.0 + denomTmp);
2304  }
2305  if (denomTmp > 0.0) {
2306  y = denomTmp;
2307  yp1 = y + 1.0;
2308  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2309  } else {
2310  sigma = 0.0;
2311  }
2312  m_dlnActCoeffMolaldT[m_indexSolvent] =
2313  2.0 /3.0 * dAdT * m_Mnaught *
2314  m_IionicMolality * sqrtI * sigma;
2315  break;
2316 
2317  case DHFORM_BETAIJ:
2318  denomTmp *= m_Aionic[0];
2319  for (size_t k = 0; k < m_kk; k++) {
2320  if (k != m_indexSolvent) {
2321  z_k = m_speciesCharge[k];
2322  m_dlnActCoeffMolaldT[k] =
2323  - z_k * z_k * numdAdTTmp / (1.0 + denomTmp);
2324  }
2325  }
2326  if (denomTmp > 0.0) {
2327  y = denomTmp;
2328  yp1 = y + 1.0;
2329  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2330  } else {
2331  sigma = 0.0;
2332  }
2333  m_dlnActCoeffMolaldT[m_indexSolvent] =
2334  2.0 /3.0 * dAdT * m_Mnaught *
2335  m_IionicMolality * sqrtI * sigma;
2336  break;
2337 
2338  case DHFORM_PITZER_BETAIJ:
2339  denomTmp *= m_Aionic[0];
2340  tmpLn = log(1.0 + denomTmp);
2341  for (size_t k = 0; k < m_kk; k++) {
2342  if (k != m_indexSolvent) {
2343  z_k = m_speciesCharge[k];
2344  m_dlnActCoeffMolaldT[k] =
2345  - z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
2346  - 2.0 * z_k * z_k * dAdT * tmpLn
2347  / (m_B_Debye * m_Aionic[0]);
2348  m_dlnActCoeffMolaldT[k] /= 3.0;
2349  }
2350  }
2351 
2352  sigma = 1.0 / (1.0 + denomTmp);
2353  m_dlnActCoeffMolaldT[m_indexSolvent] =
2354  2.0 /3.0 * dAdT * m_Mnaught *
2355  m_IionicMolality * sqrtI * sigma;
2356  break;
2357 
2358  default:
2359  printf("ERROR\n");
2360  exit(EXIT_FAILURE);
2361  break;
2362  }
2363 
2364 
2365 }
2366 
2367 /*
2368  * s_update_d2lnMolalityActCoeff_dT2() (private, const )
2369  *
2370  * Using internally stored values, this function calculates
2371  * the temperature 2nd derivative of the logarithm of the
2372  * activity coefficient
2373  * for all species in the mechanism.
2374  *
2375  * We assume that the activity coefficients are current.
2376  *
2377  * solvent activity coefficient is on the molality
2378  * scale. Its derivatives are too.
2379  */
2380 void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
2381 {
2382  double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
2383  double dAdT = dA_DebyedT_TP();
2384  double d2AdT2 = d2A_DebyedT2_TP();
2385  if (d2AdT2 == 0.0 && dAdT == 0.0) {
2386  for (size_t k = 0; k < m_kk; k++) {
2387  m_d2lnActCoeffMolaldT2[k] = 0.0;
2388  }
2389  return;
2390  }
2391 
2392  /*
2393  * Calculate a safe value for the mole fraction
2394  * of the solvent
2395  */
2396  double xmolSolvent = moleFraction(m_indexSolvent);
2397  xmolSolvent = std::max(8.689E-3, xmolSolvent);
2398 
2399 
2400  double sqrtI = sqrt(m_IionicMolality);
2401  double numd2AdT2Tmp = d2AdT2 * sqrtI;
2402  double denomTmp = m_B_Debye * sqrtI;
2403 
2404  switch (m_formDH) {
2405  case DHFORM_DILUTE_LIMIT:
2406  for (size_t k = 0; k < m_kk; k++) {
2407  m_d2lnActCoeffMolaldT2[k] =
2408  m_lnActCoeffMolal[k] * d2AdT2 / m_A_Debye;
2409  }
2410  break;
2411 
2412  case DHFORM_BDOT_AK:
2413  for (size_t k = 0; k < m_kk; k++) {
2414  z_k = m_speciesCharge[k];
2415  m_d2lnActCoeffMolaldT2[k] =
2416  - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp * m_Aionic[k]);
2417  }
2418 
2419  m_d2lnActCoeffMolaldT2[m_indexSolvent] = 0.0;
2420 
2421  coeff = 2.0 / 3.0 * d2AdT2 * m_Mnaught * sqrtI;
2422  tmp = 0.0;
2423  if (denomTmp > 0.0) {
2424  for (size_t k = 0; k < m_kk; k++) {
2425  y = denomTmp * m_Aionic[k];
2426  yp1 = y + 1.0;
2427  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2428  z_k = m_speciesCharge[k];
2429  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
2430  }
2431  }
2432  m_d2lnActCoeffMolaldT2[m_indexSolvent] += coeff * tmp;
2433  break;
2434 
2435  case DHFORM_BDOT_ACOMMON:
2436  denomTmp *= m_Aionic[0];
2437  for (size_t k = 0; k < m_kk; k++) {
2438  z_k = m_speciesCharge[k];
2439  m_d2lnActCoeffMolaldT2[k] =
2440  - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp);
2441  }
2442  if (denomTmp > 0.0) {
2443  y = denomTmp;
2444  yp1 = y + 1.0;
2445  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2446  } else {
2447  sigma = 0.0;
2448  }
2449  m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2450  2.0 /3.0 * d2AdT2 * m_Mnaught *
2451  m_IionicMolality * sqrtI * sigma;
2452  break;
2453 
2454  case DHFORM_BETAIJ:
2455  denomTmp *= m_Aionic[0];
2456  for (size_t k = 0; k < m_kk; k++) {
2457  if (k != m_indexSolvent) {
2458  z_k = m_speciesCharge[k];
2459  m_d2lnActCoeffMolaldT2[k] =
2460  - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp);
2461  }
2462  }
2463  if (denomTmp > 0.0) {
2464  y = denomTmp;
2465  yp1 = y + 1.0;
2466  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 -2.0*log(yp1));
2467  } else {
2468  sigma = 0.0;
2469  }
2470  m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2471  2.0 /3.0 * d2AdT2 * m_Mnaught *
2472  m_IionicMolality * sqrtI * sigma;
2473  break;
2474 
2475  case DHFORM_PITZER_BETAIJ:
2476  denomTmp *= m_Aionic[0];
2477  tmpLn = log(1.0 + denomTmp);
2478  for (size_t k = 0; k < m_kk; k++) {
2479  if (k != m_indexSolvent) {
2480  z_k = m_speciesCharge[k];
2481  m_d2lnActCoeffMolaldT2[k] =
2482  - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
2483  - 2.0 * z_k * z_k * d2AdT2 * tmpLn
2484  / (m_B_Debye * m_Aionic[0]);
2485  m_d2lnActCoeffMolaldT2[k] /= 3.0;
2486  }
2487  }
2488 
2489  sigma = 1.0 / (1.0 + denomTmp);
2490  m_d2lnActCoeffMolaldT2[m_indexSolvent] =
2491  2.0 /3.0 * d2AdT2 * m_Mnaught *
2492  m_IionicMolality * sqrtI * sigma;
2493  break;
2494 
2495  default:
2496  printf("ERROR\n");
2497  exit(EXIT_FAILURE);
2498  break;
2499  }
2500 }
2501 
2502 /*
2503  * s_update_dlnMolalityActCoeff_dP() (private, const )
2504  *
2505  * Using internally stored values, this function calculates
2506  * the pressure derivative of the logarithm of the
2507  * activity coefficient for all species in the mechanism.
2508  *
2509  * We assume that the activity coefficients, molalities,
2510  * and A_Debye are current.
2511  *
2512  * solvent activity coefficient is on the molality
2513  * scale. Its derivatives are too.
2514  */
2515 void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
2516 {
2517  double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
2518  int est;
2519  double dAdP = dA_DebyedP_TP();
2520  if (dAdP == 0.0) {
2521  for (size_t k = 0; k < m_kk; k++) {
2522  m_dlnActCoeffMolaldP[k] = 0.0;
2523  }
2524  return;
2525  }
2526  /*
2527  * Calculate a safe value for the mole fraction
2528  * of the solvent
2529  */
2530  double xmolSolvent = moleFraction(m_indexSolvent);
2531  xmolSolvent = std::max(8.689E-3, xmolSolvent);
2532 
2533 
2534  double sqrtI = sqrt(m_IionicMolality);
2535  double numdAdPTmp = dAdP * sqrtI;
2536  double denomTmp = m_B_Debye * sqrtI;
2537 
2538  switch (m_formDH) {
2539  case DHFORM_DILUTE_LIMIT:
2540  for (size_t k = 0; k < m_kk; k++) {
2541  m_dlnActCoeffMolaldP[k] =
2542  m_lnActCoeffMolal[k] * dAdP / m_A_Debye;
2543  }
2544  break;
2545 
2546  case DHFORM_BDOT_AK:
2547  for (size_t k = 0; k < m_kk; k++) {
2548  est = m_electrolyteSpeciesType[k];
2549  if (est == cEST_nonpolarNeutral) {
2550  m_lnActCoeffMolal[k] = 0.0;
2551  } else {
2552  z_k = m_speciesCharge[k];
2553  m_dlnActCoeffMolaldP[k] =
2554  - z_k * z_k * numdAdPTmp / (1.0 + denomTmp * m_Aionic[k]);
2555  }
2556  }
2557 
2558  m_dlnActCoeffMolaldP[m_indexSolvent] = 0.0;
2559 
2560  coeff = 2.0 / 3.0 * dAdP * m_Mnaught * sqrtI;
2561  tmp = 0.0;
2562  if (denomTmp > 0.0) {
2563  for (size_t k = 0; k < m_kk; k++) {
2564  y = denomTmp * m_Aionic[k];
2565  yp1 = y + 1.0;
2566  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2567  z_k = m_speciesCharge[k];
2568  tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
2569  }
2570  }
2571  m_dlnActCoeffMolaldP[m_indexSolvent] += coeff * tmp;
2572  break;
2573 
2574  case DHFORM_BDOT_ACOMMON:
2575  denomTmp *= m_Aionic[0];
2576  for (size_t k = 0; k < m_kk; k++) {
2577  z_k = m_speciesCharge[k];
2578  m_dlnActCoeffMolaldP[k] =
2579  - z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
2580  }
2581  if (denomTmp > 0.0) {
2582  y = denomTmp;
2583  yp1 = y + 1.0;
2584  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2585  } else {
2586  sigma = 0.0;
2587  }
2588  m_dlnActCoeffMolaldP[m_indexSolvent] =
2589  2.0 /3.0 * dAdP * m_Mnaught *
2590  m_IionicMolality * sqrtI * sigma;
2591  break;
2592 
2593  case DHFORM_BETAIJ:
2594  denomTmp *= m_Aionic[0];
2595  for (size_t k = 0; k < m_kk; k++) {
2596  if (k != m_indexSolvent) {
2597  z_k = m_speciesCharge[k];
2598  m_dlnActCoeffMolaldP[k] =
2599  - z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
2600  }
2601  }
2602  if (denomTmp > 0.0) {
2603  y = denomTmp;
2604  yp1 = y + 1.0;
2605  sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
2606  } else {
2607  sigma = 0.0;
2608  }
2609  m_dlnActCoeffMolaldP[m_indexSolvent] =
2610  2.0 /3.0 * dAdP * m_Mnaught *
2611  m_IionicMolality * sqrtI * sigma;
2612  break;
2613 
2614  case DHFORM_PITZER_BETAIJ:
2615  denomTmp *= m_Aionic[0];
2616  tmpLn = log(1.0 + denomTmp);
2617  for (size_t k = 0; k < m_kk; k++) {
2618  if (k != m_indexSolvent) {
2619  z_k = m_speciesCharge[k];
2620  m_dlnActCoeffMolaldP[k] =
2621  - z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
2622  - 2.0 * z_k * z_k * dAdP * tmpLn
2623  / (m_B_Debye * m_Aionic[0]);
2624  m_dlnActCoeffMolaldP[k] /= 3.0;
2625  }
2626  }
2627 
2628  sigma = 1.0 / (1.0 + denomTmp);
2629  m_dlnActCoeffMolaldP[m_indexSolvent] =
2630  2.0 /3.0 * dAdP * m_Mnaught *
2631  m_IionicMolality * sqrtI * sigma;
2632  break;
2633 
2634  default:
2635  printf("ERROR\n");
2636  exit(EXIT_FAILURE);
2637  break;
2638  }
2639 }
2640 
2641 /*
2642  * Updates the standard state thermodynamic functions at the current T and P of the solution.
2643  *
2644  * @internal
2645  *
2646  * This function gets called for every call to functions in this
2647  * class. It checks to see whether the temperature or pressure has changed and
2648  * thus the ss thermodynamics functions for all of the species
2649  * must be recalculated.
2650  */
2651 // void DebyeHuckel::_updateStandardStateThermo() const {
2652 // doublereal tnow = temperature();
2653 // doublereal pnow = m_Pcurrent;
2654 // if (m_waterSS) {
2655 // m_waterSS->setTempPressure(tnow, pnow);
2656 // }
2657 // m_VPSS_ptr->setState_TP(tnow, pnow);
2658 // VPStandardStateTP::updateStandardStateThermo();
2659 
2660 //}
2661 
2662 }
Generated by   doxygen 1.8.2