Cantera  2.4.0
HMWSoln.h
Go to the documentation of this file.
1 /**
2  * @file HMWSoln.h
3  * Headers for the HMWSoln ThermoPhase object, which models concentrated
4  * electrolyte solutions
5  * (see \ref thermoprops and \link Cantera::HMWSoln HMWSoln \endlink) .
6  *
7  * Class HMWSoln represents a concentrated liquid electrolyte phase which
8  * obeys the Pitzer formulation for nonideality using molality-based
9  * standard states.
10  */
11 
12 // This file is part of Cantera. See License.txt in the top-level directory or
13 // at http://www.cantera.org/license.txt for license and copyright information.
14 
15 #ifndef CT_HMWSOLN_H
16 #define CT_HMWSOLN_H
17 
18 #include "MolalityVPSSTP.h"
19 #include "cantera/base/Array.h"
20 
21 namespace Cantera
22 {
23 
24 /*!
25  * @name Temperature Dependence of the Pitzer Coefficients
26  *
27  * Note, the temperature dependence of the Gibbs free energy also depends on the
28  * temperature dependence of the standard state and the temperature dependence
29  * of the Debye-Huckel constant, which includes the dielectric constant and the
30  * density. Therefore, this expression defines only part of the temperature
31  * dependence for the mixture thermodynamic functions.
32  *
33  * PITZER_TEMP_CONSTANT
34  * All coefficients are considered constant wrt temperature
35  * PITZER_TEMP_LINEAR
36  * All coefficients are assumed to have a linear dependence
37  * wrt to temperature.
38  * PITZER_TEMP_COMPLEX1
39  * All coefficients are assumed to have a complex functional
40  * based dependence wrt temperature; See:
41  * (Silvester, Pitzer, J. Phys. Chem. 81, 19 1822 (1977)).
42  *
43  * beta0 = q0 + q3(1/T - 1/Tr) + q4(ln(T/Tr)) +
44  * q1(T - Tr) + q2(T**2 - Tr**2)
45  */
46 //@{
47 #define PITZER_TEMP_CONSTANT 0
48 #define PITZER_TEMP_LINEAR 1
49 #define PITZER_TEMP_COMPLEX1 2
50 //@}
51 
52 /*
53  * @name ways to calculate the value of A_Debye
54  *
55  * These defines determine the way A_Debye is calculated
56  */
57 //@{
58 #define A_DEBYE_CONST 0
59 #define A_DEBYE_WATER 1
60 //@}
61 
62 class WaterProps;
63 
64 /**
65  * Class HMWSoln represents a dilute or concentrated liquid electrolyte
66  * phase which obeys the Pitzer formulation for nonideality.
67  *
68  * As a prerequisite to the specification of thermodynamic quantities,
69  * The concentrations of the ionic species are assumed to obey the
70  * electroneutrality condition.
71  *
72  * ## Specification of Species Standard State Properties
73  *
74  * The solvent is assumed to be liquid water. A real model for liquid water
75  * (IAPWS 1995 formulation) is used as its standard state. All standard state
76  * properties for the solvent are based on this real model for water, and
77  * involve function calls to the object that handles the real water model,
78  * #Cantera::WaterPropsIAPWS.
79  *
80  * The standard states for solutes are on the unit molality basis. Therefore, in
81  * the documentation below, the normal \f$ o \f$ superscript is replaced with
82  * the \f$ \triangle \f$ symbol. The reference state symbol is now
83  * \f$ \triangle, ref \f$.
84  *
85  * It is assumed that the reference state thermodynamics may be obtained by a
86  * pointer to a populated species thermodynamic property manager class (see
87  * ThermoPhase::m_spthermo). How to relate pressure changes to the reference
88  * state thermodynamics is resolved at this level.
89  *
90  * For solutes that rely on ThermoPhase::m_spthermo, are assumed to have an
91  * incompressible standard state mechanical property. In other words, the molar
92  * volumes are independent of temperature and pressure.
93  *
94  * For these incompressible, standard states, the molar internal energy is
95  * independent of pressure. Since the thermodynamic properties are specified by
96  * giving the standard-state enthalpy, the term \f$ P_0 \hat v\f$ is subtracted
97  * from the specified molar enthalpy to compute the molar internal energy. The
98  * entropy is assumed to be independent of the pressure.
99  *
100  * The enthalpy function is given by the following relation.
101  *
102  * \f[
103  * h^\triangle_k(T,P) = h^{\triangle,ref}_k(T)
104  * + \tilde{v}_k \left( P - P_{ref} \right)
105  * \f]
106  *
107  * For an incompressible, stoichiometric substance, the molar internal energy is
108  * independent of pressure. Since the thermodynamic properties are specified by
109  * giving the standard-state enthalpy, the term \f$ P_{ref} \tilde v\f$ is
110  * subtracted from the specified reference molar enthalpy to compute the molar
111  * internal energy.
112  *
113  * \f[
114  * u^\triangle_k(T,P) = h^{\triangle,ref}_k(T) - P_{ref} \tilde{v}_k
115  * \f]
116  *
117  * The solute standard state heat capacity and entropy are independent of
118  * pressure. The solute standard state Gibbs free energy is obtained from the
119  * enthalpy and entropy functions.
120  *
121  * The current model assumes that an incompressible molar volume for all
122  * solutes. The molar volume for the water solvent, however, is obtained from a
123  * pure water equation of state, waterSS. Therefore, the water standard state
124  * varies with both T and P. It is an error to request standard state water
125  * properties at a T and P where the water phase is not a stable phase, i.e.,
126  * beyond its spinodal curve.
127  *
128  * ## Specification of Solution Thermodynamic Properties
129  *
130  * Chemical potentials of the solutes, \f$ \mu_k \f$, and the solvent, \f$ \mu_o
131  * \f$, which are based on the molality form, have the following general format:
132  *
133  * \f[
134  * \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} \frac{m_k}{m^\triangle})
135  * \f]
136  * \f[
137  * \mu_o = \mu^o_o(T,P) + RT ln(a_o)
138  * \f]
139  *
140  * where \f$ \gamma_k^{\triangle} \f$ is the molality based activity coefficient
141  * for species \f$k\f$.
142  *
143  * Individual activity coefficients of ions can not be independently measured.
144  * Instead, only binary pairs forming electroneutral solutions can be measured.
145  * This problem leads to a redundancy in the evaluation of species standard
146  * state properties. The redundancy issue is resolved by setting the standard
147  * state chemical potential enthalpy, entropy, and volume for the hydrogen ion,
148  * H+, to zero, for every temperature and pressure. After this convention is
149  * applied, all other standard state properties of ionic species contain
150  * meaningful information.
151  *
152  * ### Ionic Strength
153  *
154  * Most of the parameterizations within the model use the ionic strength as a
155  * key variable. The ionic strength, \f$ I\f$ is defined as follows
156  *
157  * \f[
158  * I = \frac{1}{2} \sum_k{m_k z_k^2}
159  * \f]
160  *
161  * \f$ m_k \f$ is the molality of the kth species. \f$ z_k \f$ is the charge of
162  * the kth species. Note, the ionic strength is a defined units quantity. The
163  * molality has defined units of gmol kg-1, and therefore the ionic strength has
164  * units of sqrt(gmol/kg).
165  *
166  * ### Specification of the Excess Gibbs Free Energy
167  *
168  * Pitzer's formulation may best be represented as a specification of the excess
169  * Gibbs free energy, \f$ G^{ex} \f$, defined as the deviation of the total
170  * Gibbs free energy from that of an ideal molal solution.
171  * \f[
172  * G = G^{id} + G^{ex}
173  * \f]
174  *
175  * The ideal molal solution contribution, not equal to an ideal solution
176  * contribution and in fact containing a singularity at the zero solvent mole
177  * fraction limit, is given below.
178  * \f[
179  * G^{id} = n_o \mu^o_o + \sum_{k\ne o} n_k \mu_k^{\triangle}
180  * + \tilde{M}_o n_o ( RT (\sum{m_i(\ln(m_i)-1)}))
181  * \f]
182  *
183  * From the excess Gibbs free energy formulation, the activity coefficient
184  * expression and the osmotic coefficient expression for the solvent may be
185  * defined, by taking the appropriate derivatives. Using this approach
186  * guarantees that the entire system will obey the Gibbs-Duhem relations.
187  *
188  * Pitzer employs the following general expression for the excess Gibbs free
189  * energy
190  *
191  * \f[
192  * \begin{array}{cclc}
193  * \frac{G^{ex}}{\tilde{M}_o n_o RT} &= &
194  * \left( \frac{4A_{Debye}I}{3b} \right) \ln(1 + b \sqrt{I})
195  * + 2 \sum_c \sum_a m_c m_a B_{ca}
196  * + \sum_c \sum_a m_c m_a Z C_{ca}
197  * \\&&
198  * + \sum_{c < c'} \sum m_c m_{c'} \left[ 2 \Phi_{c{c'}} + \sum_a m_a \Psi_{c{c'}a} \right]
199  * + \sum_{a < a'} \sum m_a m_{a'} \left[ 2 \Phi_{a{a'}} + \sum_c m_c \Psi_{a{a'}c} \right]
200  * \\&&
201  * + 2 \sum_n \sum_c m_n m_c \lambda_{nc} + 2 \sum_n \sum_a m_n m_a \lambda_{na}
202  * + 2 \sum_{n < n'} \sum m_n m_{n'} \lambda_{n{n'}}
203  * + \sum_n m^2_n \lambda_{nn}
204  * \end{array}
205  * \f]
206  *
207  * *a* is a subscript over all anions, *c* is a subscript extending over all
208  * cations, and *i* is a subscript that extends over all anions and cations.
209  * *n* is a subscript that extends only over neutral solute molecules. The
210  * second line contains cross terms where cations affect cations and/or
211  * cation/anion pairs, and anions affect anions or cation/anion pairs. Note part
212  * of the coefficients, \f$ \Phi_{c{c'}} \f$ and \f$ \Phi_{a{a'}} \f$ stem from
213  * the theory of unsymmetrical mixing of electrolytes with different charges.
214  * This theory depends on the total ionic strength of the solution, and
215  * therefore, \f$ \Phi_{c{c'}} \f$ and \f$ \Phi_{a{a'}} \f$ will depend on
216  * *I*, the ionic strength. \f$ B_{ca}\f$ is a strong function of the
217  * total ionic strength, *I*, of the electrolyte. The rest of the coefficients
218  * are assumed to be independent of the molalities or ionic strengths. However,
219  * all coefficients are potentially functions of the temperature and pressure
220  * of the solution.
221  *
222  * *A* is the Debye-Huckel constant. Its specification is described in its
223  * own section below.
224  *
225  * \f$ I\f$ is the ionic strength of the solution, and is given by:
226  *
227  * \f[
228  * I = \frac{1}{2} \sum_k{m_k z_k^2}
229  * \f]
230  *
231  * In contrast to several other Debye-Huckel implementations (see \ref
232  * DebyeHuckel), the parameter \f$ b\f$ in the above equation is a constant that
233  * does not vary with respect to ion identity. This is an important
234  * simplification as it avoids troubles with satisfaction of the Gibbs-Duhem
235  * analysis.
236  *
237  * The function \f$ Z \f$ is given by
238  *
239  * \f[
240  * Z = \sum_i m_i \left| z_i \right|
241  * \f]
242  *
243  * The value of \f$ B_{ca}\f$ is given by the following function
244  *
245  * \f[
246  * B_{ca} = \beta^{(0)}_{ca} + \beta^{(1)}_{ca} g(\alpha^{(1)}_{ca} \sqrt{I})
247  * + \beta^{(2)}_{ca} g(\alpha^{(2)}_{ca} \sqrt{I})
248  * \f]
249  *
250  * where
251  *
252  * \f[
253  * g(x) = 2 \frac{(1 - (1 + x)\exp[-x])}{x^2}
254  * \f]
255  *
256  * The formulation for \f$ B_{ca}\f$ combined with the formulation of the Debye-
257  * Huckel term in the eqn. for the excess Gibbs free energy stems essentially
258  * from an empirical fit to the ionic strength dependent data based over a wide
259  * sampling of binary electrolyte systems. \f$ C_{ca} \f$, \f$ \lambda_{nc} \f$,
260  * \f$ \lambda_{na} \f$, \f$ \lambda_{nn} \f$, \f$ \Psi_{c{c'}a} \f$, \f$
261  * \Psi_{a{a'}c} \f$ are experimentally derived coefficients that may have
262  * pressure and/or temperature dependencies.
263  *
264  * The \f$ \Phi_{c{c'}} \f$ and \f$ \Phi_{a{a'}} \f$ formulations are slightly
265  * more complicated. \f$ b \f$ is a universal constant defined to be equal to
266  * \f$ 1.2\ kg^{1/2}\ gmol^{-1/2} \f$. The exponential coefficient \f$
267  * \alpha^{(1)}_{ca} \f$ is usually fixed at \f$ \alpha^{(1)}_{ca} = 2.0\
268  * kg^{1/2} gmol^{-1/2}\f$ except for 2-2 electrolytes, while other parameters
269  * were fit to experimental data. For 2-2 electrolytes, \f$ \alpha^{(1)}_{ca} =
270  * 1.4\ kg^{1/2}\ gmol^{-1/2}\f$ is used in combination with either \f$
271  * \alpha^{(2)}_{ca} = 12\ kg^{1/2}\ gmol^{-1/2}\f$ or \f$ \alpha^{(2)}_{ca} = k
272  * A_\psi \f$, where *k* is a constant. For electrolytes other than 2-2
273  * electrolytes the \f$ \beta^{(2)}_{ca} g(\alpha^{(2)}_{ca} \sqrt{I}) \f$ term
274  * is not used in the fitting procedure; it is only used for divalent metal
275  * solfates and other high-valence electrolytes which exhibit significant
276  * association at low ionic strengths.
277  *
278  * The \f$ \beta^{(0)}_{ca} \f$, \f$ \beta^{(1)}_{ca}\f$, \f$ \beta^{(2)}_{ca}
279  * \f$, and \f$ C_{ca} \f$ binary coefficients are referred to as ion-
280  * interaction or Pitzer parameters. These Pitzer parameters may vary with
281  * temperature and pressure but they do not depend on the ionic strength. Their
282  * values and temperature derivatives of their values have been tabulated for a
283  * range of electrolytes
284  *
285  * The \f$ \Phi_{c{c'}} \f$ and \f$ \Phi_{a{a'}} \f$ contributions, which
286  * capture cation-cation and anion-anion interactions, also have an ionic
287  * strength dependence.
288  *
289  * Ternary contributions \f$ \Psi_{c{c'}a} \f$ and \f$ \Psi_{a{a'}c} \f$ have
290  * been measured also for some systems. The success of the Pitzer method lies in
291  * its ability to model nonlinear activity coefficients of complex
292  * multicomponent systems with just binary and minor ternary contributions,
293  * which can be independently measured in binary or ternary subsystems.
294  *
295  * ### Multicomponent Activity Coefficients for Solutes
296  *
297  * The formulas for activity coefficients of solutes may be obtained by taking
298  * the following derivative of the excess Gibbs Free Energy formulation
299  * described above:
300  *
301  * \f[
302  * \ln(\gamma_k^\triangle) = \frac{d\left( \frac{G^{ex}}{M_o n_o RT} \right)}{d(m_k)}\Bigg|_{n_i}
303  * \f]
304  *
305  * In the formulas below the following conventions are used. The subscript *M*
306  * refers to a particular cation. The subscript X refers to a particular anion,
307  * whose activity is being currently evaluated. the subscript *a* refers to a
308  * summation over all anions in the solution, while the subscript *c* refers to
309  * a summation over all cations in the solutions.
310  *
311  * The activity coefficient for a particular cation *M* is given by
312  *
313  * \f[
314  * \ln(\gamma_M^\triangle) = -z_M^2(F) + \sum_a m_a \left( 2 B_{Ma} + Z C_{Ma} \right)
315  * + z_M \left( \sum_a \sum_c m_a m_c C_{ca} \right)
316  * + \sum_c m_c \left[ 2 \Phi_{Mc} + \sum_a m_a \Psi_{Mca} \right]
317  * + \sum_{a < a'} \sum m_a m_{a'} \Psi_{Ma{a'}}
318  * + 2 \sum_n m_n \lambda_{nM}
319  * \f]
320  *
321  * The activity coefficient for a particular anion *X* is given by
322  *
323  * \f[
324  * \ln(\gamma_X^\triangle) = -z_X^2(F) + \sum_a m_c \left( 2 B_{cX} + Z C_{cX} \right)
325  * + \left|z_X \right| \left( \sum_a \sum_c m_a m_c C_{ca} \right)
326  * + \sum_a m_a \left[ 2 \Phi_{Xa} + \sum_c m_c \Psi_{cXa} \right]
327  * + \sum_{c < c'} \sum m_c m_{c'} \Psi_{c{c'}X}
328  * + 2 \sum_n m_n \lambda_{nM}
329  * \f]
330  * where the function \f$ F \f$ is given by
331  *
332  * \f[
333  * F = - A_{\phi} \left[ \frac{\sqrt{I}}{1 + b \sqrt{I}}
334  * + \frac{2}{b} \ln{\left(1 + b\sqrt{I}\right)} \right]
335  * + \sum_a \sum_c m_a m_c B'_{ca}
336  * + \sum_{c < c'} \sum m_c m_{c'} \Phi'_{c{c'}}
337  * + \sum_{a < a'} \sum m_a m_{a'} \Phi'_{a{a'}}
338  * \f]
339  *
340  * We have employed the definition of \f$ A_{\phi} \f$, also used by Pitzer
341  * which is equal to
342  *
343  * \f[
344  * A_{\phi} = \frac{A_{Debye}}{3}
345  * \f]
346  *
347  * In the above formulas, \f$ \Phi'_{c{c'}} \f$ and \f$ \Phi'_{a{a'}} \f$ are the
348  * ionic strength derivatives of \f$ \Phi_{c{c'}} \f$ and \f$ \Phi_{a{a'}} \f$,
349  * respectively.
350  *
351  * The function \f$ B'_{MX} \f$ is defined as:
352  *
353  * \f[
354  * B'_{MX} = \left( \frac{\beta^{(1)}_{MX} h(\alpha^{(1)}_{MX} \sqrt{I})}{I} \right)
355  * \left( \frac{\beta^{(2)}_{MX} h(\alpha^{(2)}_{MX} \sqrt{I})}{I} \right)
356  * \f]
357  *
358  * where \f$ h(x) \f$ is defined as
359  *
360  * \f[
361  * h(x) = g'(x) \frac{x}{2} =
362  * \frac{2\left(1 - \left(1 + x + \frac{x^2}{2} \right)\exp(-x) \right)}{x^2}
363  * \f]
364  *
365  * The activity coefficient for neutral species *N* is given by
366  *
367  * \f[
368  * \ln(\gamma_N^\triangle) = 2 \left( \sum_i m_i \lambda_{iN}\right)
369  * \f]
370  *
371  * ### Activity of the Water Solvent
372  *
373  * The activity for the solvent water,\f$ a_o \f$, is not independent and must
374  * be determined either from the Gibbs-Duhem relation or from taking the
375  * appropriate derivative of the same excess Gibbs free energy function as was
376  * used to formulate the solvent activity coefficients. Pitzer's description
377  * follows the later approach to derive a formula for the osmotic coefficient,
378  * \f$ \phi \f$.
379  *
380  * \f[
381  * \phi - 1 = - \left( \frac{d\left(\frac{G^{ex}}{RT} \right)}{d(\tilde{M}_o n_o)} \right)
382  * \frac{1}{\sum_{i \ne 0} m_i}
383  * \f]
384  *
385  * The osmotic coefficient may be related to the water activity by the following relation:
386  *
387  * \f[
388  * \phi = - \frac{1}{\tilde{M}_o \sum_{i \neq o} m_i} \ln(a_o)
389  * = - \frac{n_o}{\sum_{i \neq o}n_i} \ln(a_o)
390  * \f]
391  *
392  * The result is the following
393  *
394  * \f[
395  * \begin{array}{ccclc}
396  * \phi - 1 &= &
397  * \frac{2}{\sum_{i \ne 0} m_i}
398  * \bigg[ &
399  * - A_{\phi} \frac{I^{3/2}}{1 + b \sqrt{I}}
400  * + \sum_c \sum_a m_c m_a \left( B^{\phi}_{ca} + Z C_{ca}\right)
401  * \\&&&
402  * + \sum_{c < c'} \sum m_c m_{c'} \left[ \Phi^{\phi}_{c{c'}} + \sum_a m_a \Psi_{c{c'}a} \right]
403  * + \sum_{a < a'} \sum m_a m_{a'} \left[ \Phi^{\phi}_{a{a'}} + \sum_c m_c \Psi_{a{a'}c} \right]
404  * \\&&&
405  * + \sum_n \sum_c m_n m_c \lambda_{nc} + \sum_n \sum_a m_n m_a \lambda_{na}
406  * + \sum_{n < n'} \sum m_n m_{n'} \lambda_{n{n'}}
407  * + \frac{1}{2} \left( \sum_n m^2_n \lambda_{nn}\right)
408  * \bigg]
409  * \end{array}
410  * \f]
411  *
412  * It can be shown that the expression
413  *
414  * \f[
415  * B^{\phi}_{ca} = \beta^{(0)}_{ca} + \beta^{(1)}_{ca} \exp{(- \alpha^{(1)}_{ca} \sqrt{I})}
416  * + \beta^{(2)}_{ca} \exp{(- \alpha^{(2)}_{ca} \sqrt{I} )}
417  * \f]
418  *
419  * is consistent with the expression \f$ B_{ca} \f$ in the \f$ G^{ex} \f$
420  * expression after carrying out the derivative wrt \f$ m_M \f$.
421  *
422  * Also taking into account that \f$ {\Phi}_{c{c'}} \f$ and
423  * \f$ {\Phi}_{a{a'}} \f$ has an ionic strength dependence.
424  *
425  * \f[
426  * \Phi^{\phi}_{c{c'}} = {\Phi}_{c{c'}} + I \frac{d{\Phi}_{c{c'}}}{dI}
427  * \f]
428  *
429  * \f[
430  * \Phi^{\phi}_{a{a'}} = \Phi_{a{a'}} + I \frac{d\Phi_{a{a'}}}{dI}
431  * \f]
432  *
433  * ### Temperature and Pressure Dependence of the Pitzer Parameters
434  *
435  * In general most of the coefficients introduced in the previous section may
436  * have a temperature and pressure dependence. The temperature and pressure
437  * dependence of these coefficients strongly influence the value of the excess
438  * Enthalpy and excess Volumes of Pitzer solutions. Therefore, these are readily
439  * measurable quantities. HMWSoln provides several different methods for putting
440  * these dependencies into the coefficients. HMWSoln has an implementation
441  * described by Silverter and Pitzer (1977), which was used to fit experimental
442  * data for NaCl over an extensive range, below the critical temperature of
443  * water. They found a temperature functional form for fitting the 3 following
444  * coefficients that describe the Pitzer parameterization for a single salt to
445  * be adequate to describe how the excess Gibbs free energy values for the
446  * binary salt changes with respect to temperature. The following functional
447  * form was used to fit the temperature dependence of the Pitzer Coefficients
448  * for each cation - anion pair, M X.
449  *
450  * \f[
451  * \beta^{(0)}_{MX} = q^{b0}_0
452  * + q^{b0}_1 \left( T - T_r \right)
453  * + q^{b0}_2 \left( T^2 - T_r^2 \right)
454  * + q^{b0}_3 \left( \frac{1}{T} - \frac{1}{T_r}\right)
455  * + q^{b0}_4 \ln \left( \frac{T}{T_r} \right)
456  * \f]
457  * \f[
458  * \beta^{(1)}_{MX} = q^{b1}_0 + q^{b1}_1 \left( T - T_r \right)
459  * + q^{b1}_{2} \left( T^2 - T_r^2 \right)
460  * \f]
461  * \f[
462  * C^{\phi}_{MX} = q^{Cphi}_0
463  * + q^{Cphi}_1 \left( T - T_r \right)
464  * + q^{Cphi}_2 \left( T^2 - T_r^2 \right)
465  * + q^{Cphi}_3 \left( \frac{1}{T} - \frac{1}{T_r}\right)
466  * + q^{Cphi}_4 \ln \left( \frac{T}{T_r} \right)
467  * \f]
468  *
469  * where
470  *
471  * \f[
472  * C^{\phi}_{MX} = 2 {\left| z_M z_X \right|}^{1/2} C_{MX}
473  * \f]
474  *
475  * In later papers, Pitzer has added additional temperature dependencies to all
476  * of the other remaining second and third order virial coefficients. Some of
477  * these dependencies are justified and motivated by theory. Therefore, a
478  * formalism wherein all of the coefficients in the base theory have temperature
479  * dependencies associated with them has been implemented within the HMWSoln
480  * object. Much of the formalism, however, has been unexercised.
481  *
482  * In the HMWSoln object, the temperature dependence of the Pitzer parameters
483  * are specified in the following way.
484  *
485  * - PIZTER_TEMP_CONSTANT - string name "CONSTANT"
486  * - Assumes that all coefficients are independent of temperature
487  * and pressure
488  * - PIZTER_TEMP_COMPLEX1 - string name "COMPLEX" or "COMPLEX1"
489  * - Uses the full temperature dependence for the
490  * \f$\beta^{(0)}_{MX} \f$ (5 coeffs),
491  * the \f$\beta^{(1)}_{MX} \f$ (3 coeffs),
492  * and \f$ C^{\phi}_{MX} \f$ (5 coeffs) parameters described above.
493  * - PITZER_TEMP_LINEAR - string name "LINEAR"
494  * - Uses just the temperature dependence for the
495  * \f$\beta^{(0)}_{MX} \f$, the \f$\beta^{(1)}_{MX} \f$,
496  * and \f$ C^{\phi}_{MX} \f$ coefficients described above.
497  * There are 2 coefficients for each term.
498  *
499  * The temperature dependence is specified in an attributes field in the
500  * `activityCoefficients` XML block, called `TempModel`. Permissible values for
501  * that attribute are `CONSTANT`, `COMPLEX1`, and `LINEAR`.
502  *
503  * The specification of the binary interaction between a cation and an anion is
504  * given by the coefficients, \f$ B_{MX}\f$ and \f$ C_{MX}\f$ The specification
505  * of \f$ B_{MX}\f$ is a function of \f$\beta^{(0)}_{MX} \f$,
506  * \f$\beta^{(1)}_{MX} \f$, \f$\beta^{(2)}_{MX} \f$, \f$\alpha^{(1)}_{MX} \f$,
507  * and \f$\alpha^{(2)}_{MX} \f$. \f$ C_{MX}\f$ is calculated from
508  * \f$C^{\phi}_{MX} \f$ from the formula above. All of the underlying
509  * coefficients are specified in the XML element block `binarySaltParameters`,
510  * which has the attribute `cation` and `anion` to identify the interaction. XML
511  * elements named `beta0, beta1, beta2, Cphi, Alpha1, Alpha2` within each
512  * `binarySaltParameters` block specify the parameters. Within each of these
513  * blocks multiple parameters describing temperature or pressure dependence are
514  * serially listed in the order that they appear in the equation in this
515  * document. An example of the `beta0` block that fits the `COMPLEX1`
516  * temperature dependence given above is
517  *
518  * @code
519  * <binarySaltParameters cation="Na+" anion="OH-">
520  * <beta0> q0, q1, q2, q3, q4 </beta0>
521  * </binarySaltParameters>
522  * @endcode
523  *
524  * The parameters for \f$ \beta^{(0)}\f$ fit the following equation:
525  *
526  * \f[
527  * \beta^{(0)} = q_0^{{\beta}0} + q_1^{{\beta}0} \left( T - T_r \right)
528  * + q_2^{{\beta}0} \left( T^2 - T_r^2 \right)
529  * + q_3^{{\beta}0} \left( \frac{1}{T} - \frac{1}{T_r} \right)
530  * + q_4^{{\beta}0} \ln \left( \frac{T}{T_r} \right)
531  * \f]
532  *
533  * This same `COMPLEX1` temperature dependence given above is used for the
534  * following parameters:
535  * \f$ \beta^{(0)}_{MX} \f$, \f$ \beta^{(1)}_{MX} \f$,
536  * \f$ \beta^{(2)}_{MX} \f$, \f$ \Theta_{cc'} \f$, \f$\Theta_{aa'} \f$,
537  * \f$ \Psi_{c{c'}a} \f$ and \f$ \Psi_{ca{a'}} \f$.
538  *
539  * ### Like-Charged Binary Ion Parameters and the Mixing Parameters
540  *
541  * The previous section contained the functions, \f$ \Phi_{c{c'}} \f$,
542  * \f$ \Phi_{a{a'}} \f$ and their derivatives wrt the ionic strength, \f$
543  * \Phi'_{c{c'}} \f$ and \f$ \Phi'_{a{a'}} \f$. Part of these terms come from
544  * theory.
545  *
546  * Since like charged ions repel each other and are generally not near each
547  * other, the virial coefficients for same-charged ions are small. However,
548  * Pitzer doesn't ignore these in his formulation. Relatively larger and longer
549  * range terms between like-charged ions exist however, which appear only for
550  * unsymmetrical mixing of same-sign charged ions with different charges. \f$
551  * \Phi_{ij} \f$, where \f$ ij \f$ is either \f$ a{a'} \f$ or \f$ c{c'} \f$ is
552  * given by
553  *
554  * \f[
555  * {\Phi}_{ij} = \Theta_{ij} + \,^E \Theta_{ij}(I)
556  * \f]
557  *
558  * \f$ \Theta_{ij} \f$ is the small virial coefficient expansion term. Dependent
559  * in general on temperature and pressure, its ionic strength dependence is
560  * ignored in Pitzer's approach. \f$ \,^E\Theta_{ij}(I) \f$ accounts for the
561  * electrostatic unsymmetrical mixing effects and is dependent only on the
562  * charges of the ions i, j, the total ionic strength and on the dielectric
563  * constant and density of the solvent. This seems to be a relatively well-
564  * documented part of the theory. They theory below comes from Pitzer summation
565  * (Pitzer) in the appendix. It's also mentioned in Bethke's book (Bethke), and
566  * the equations are summarized in Harvie & Weare (1980). Within the code, \f$
567  * \,^E\Theta_{ij}(I) \f$ is evaluated according to the algorithm described in
568  * Appendix B [Pitzer] as
569  *
570  * \f[
571  * \,^E\Theta_{ij}(I) = \left( \frac{z_i z_j}{4I} \right)
572  * \left( J(x_{ij}) - \frac{1}{2} J(x_{ii})
573  * - \frac{1}{2} J(x_{jj}) \right)
574  * \f]
575  *
576  * where \f$ x_{ij} = 6 z_i z_j A_{\phi} \sqrt{I} \f$ and
577  *
578  * \f[
579  * J(x) = \frac{1}{x} \int_0^{\infty}{\left( 1 + q +
580  * \frac{1}{2} q^2 - e^q \right) y^2 dy}
581  * \f]
582  *
583  * and \f$ q = - (\frac{x}{y}) e^{-y} \f$. \f$ J(x) \f$ is evaluated by
584  * numerical integration.
585  *
586  * The \f$ \Theta_{ij} \f$ term is a constant that is specified by the XML
587  * element `thetaCation` and `thetaAnion`, which has the attribute `cation1`,
588  * `cation2` and `anion1`, `anion2` respectively to identify the interaction. No
589  * temperature or pressure dependence of this parameter is currently allowed. An
590  * example of the block is presented below.
591  *
592  * @code
593  * <thetaCation cation1="Na+" cation2="H+">
594  * <Theta> 0.036 </Theta>
595  * </thetaCation>
596  * @endcode
597  *
598  * ### Ternary Pitzer Parameters
599  *
600  * The \f$ \Psi_{c{c'}a} \f$ and \f$ \Psi_{ca{a'}} \f$ terms represent ternary
601  * interactions between two cations and an anion and two anions and a cation,
602  * respectively. In Pitzer's implementation these terms are usually small in
603  * absolute size. Currently these parameters do not have any dependence on
604  * temperature, pressure, or ionic strength.
605  *
606  * Their values are input using the XML element `psiCommonCation` and
607  * `psiCommonAnion`. The species id's are specified in attribute fields in the
608  * XML element. The fields `cation`, `anion1`, and `anion2` are used for
609  * `psiCommonCation`. The fields `anion`, `cation1` and `cation2` are used for
610  * `psiCommonAnion`. An example block is given below. The `Theta` field below is
611  * a duplicate of the `thetaAnion` field mentioned above. The two fields are
612  * input into the same block for convenience, and because their data are highly
613  * correlated, in practice. It is an error for the two blocks to specify
614  * different information about thetaAnion (or thetaCation) in different blocks.
615  * It's ok to specify duplicate but consistent information in multiple blocks.
616  *
617  * @code
618  * <psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
619  * <Theta> -0.05 </Theta>
620  * <Psi> -0.006 </Psi>
621  * </psiCommonCation>
622  * @endcode
623  *
624  * ### Treatment of Neutral Species
625  *
626  * Binary virial-coefficient-like interactions between two neutral species may
627  * be specified in the \f$ \lambda_{mn} \f$ terms that appear in the formulas
628  * above. Currently these interactions are independent of temperature, pressure,
629  * and ionic strength. Also, currently, the neutrality of the species are not
630  * checked. Therefore, this interaction may involve charged species in the
631  * solution as well. The identity of the species is specified by the `species1`
632  * and `species2` attributes to the XML `lambdaNeutral` node. These terms are
633  * symmetrical; `species1` and `species2` may be reversed and the term will be
634  * the same. An example is given below.
635  *
636  * @code
637  * <lambdaNeutral species1="CO2" species2="CH4">
638  * <lambda> 0.05 </lambda>
639  * </lambdaNeutral>
640  * @endcode
641  *
642  * ## Example of the Specification of Parameters for the Activity Coefficients
643  *
644  * An example is given below.
645  *
646  * An example `activityCoefficients` XML block for this formulation is supplied
647  * below
648  *
649  * @code
650  * <activityCoefficients model="Pitzer" TempModel="complex1">
651  * <!-- Pitzer Coefficients
652  * These coefficients are from Pitzer's main
653  * paper, in his book.
654  * -->
655  * <A_Debye model="water" />
656  * <ionicRadius default="3.042843" units="Angstroms">
657  * </ionicRadius>
658  * <binarySaltParameters cation="Na+" anion="Cl-">
659  * <beta0> 0.0765, 0.008946, -3.3158E-6,
660  * -777.03, -4.4706
661  * </beta0>
662  * <beta1> 0.2664, 6.1608E-5, 1.0715E-6, 0.0, 0.0 </beta1>
663  * <beta2> 0.0, 0.0, 0.0, 0.0, 0.0 </beta2>
664  * <Cphi> 0.00127, -4.655E-5, 0.0,
665  * 33.317, 0.09421
666  * </Cphi>
667  * <Alpha1> 2.0 </Alpha1>
668  * </binarySaltParameters>
669  *
670  * <binarySaltParameters cation="H+" anion="Cl-">
671  * <beta0> 0.1775, 0.0, 0.0, 0.0, 0.0 </beta0>
672  * <beta1> 0.2945, 0.0, 0.0, 0.0, 0.0 </beta1>
673  * <beta2> 0.0, 0.0, 0.0, 0.0, 0.0 </beta2>
674  * <Cphi> 0.0008, 0.0, 0.0, 0.0, 0.0 </Cphi>
675  * <Alpha1> 2.0 </Alpha1>
676  * </binarySaltParameters>
677  *
678  * <binarySaltParameters cation="Na+" anion="OH-">
679  * <beta0> 0.0864, 0.0, 0.0, 0.0, 0.0 </beta0>
680  * <beta1> 0.253, 0.0, 0.0 0.0, 0.0 </beta1>
681  * <beta2> 0.0 0.0, 0.0, 0.0, 0.0 </beta2>
682  * <Cphi> 0.0044, 0.0, 0.0, 0.0, 0.0 </Cphi>
683  * <Alpha1> 2.0 </Alpha1>
684  * </binarySaltParameters>
685  *
686  * <thetaAnion anion1="Cl-" anion2="OH-">
687  * <Theta> -0.05, 0.0, 0.0, 0.0, 0.0 </Theta>
688  * </thetaAnion>
689  *
690  * <psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
691  * <Theta> -0.05, 0.0, 0.0, 0.0, 0.0 </Theta>
692  * <Psi> -0.006 </Psi>
693  * </psiCommonCation>
694  *
695  * <thetaCation cation1="Na+" cation2="H+">
696  * <Theta> 0.036, 0.0, 0.0, 0.0, 0.0 </Theta>
697  * </thetaCation>
698  *
699  * <psiCommonAnion anion="Cl-" cation1="Na+" cation2="H+">
700  * <Theta> 0.036, 0.0, 0.0, 0.0, 0.0 </Theta>
701  * <Psi> -0.004 </Psi>
702  * </psiCommonAnion>
703  * </activityCoefficients>
704  * @endcode
705  *
706  * ### Specification of the Debye-Huckel Constant
707  *
708  * In the equations above, the formula for \f$ A_{Debye} \f$ is needed. The
709  * HMWSoln object uses two methods for specifying these quantities. The default
710  * method is to assume that \f$ A_{Debye} \f$ is a constant, given in the
711  * initialization process, and stored in the member double, m_A_Debye.
712  * Optionally, a full water treatment may be employed that makes
713  * \f$ A_{Debye} \f$ a full function of *T* and *P* and creates nontrivial
714  * entries for the excess heat capacity, enthalpy, and excess volumes of
715  * solution.
716  *
717  * \f[
718  * A_{Debye} = \frac{F e B_{Debye}}{8 \pi \epsilon R T} {\left( C_o \tilde{M}_o \right)}^{1/2}
719  * \f]
720  * where
721  *
722  * \f[
723  * B_{Debye} = \frac{F} {{(\frac{\epsilon R T}{2})}^{1/2}}
724  * \f]
725  * Therefore:
726  * \f[
727  * A_{Debye} = \frac{1}{8 \pi}
728  * {\left(\frac{2 N_a \rho_o}{1000}\right)}^{1/2}
729  * {\left(\frac{N_a e^2}{\epsilon R T }\right)}^{3/2}
730  * \f]
731  *
732  * Units = sqrt(kg/gmol)
733  *
734  * where
735  * - \f$ N_a \f$ is Avogadro's number
736  * - \f$ \rho_w \f$ is the density of water
737  * - \f$ e \f$ is the electronic charge
738  * - \f$ \epsilon = K \epsilon_o \f$ is the permittivity of water
739  * - \f$ K \f$ is the dielectric constant of water,
740  * - \f$ \epsilon_o \f$ is the permittivity of free space.
741  * - \f$ \rho_o \f$ is the density of the solvent in its standard state.
742  *
743  * Nominal value at 298 K and 1 atm = 1.172576 (kg/gmol)^(1/2)
744  * based on:
745  * - \f$ \epsilon / \epsilon_0 \f$ = 78.54 (water at 25C)
746  * - T = 298.15 K
747  * - B_Debye = 3.28640E9 (kg/gmol)^(1/2) / m
748  *
749  * An example of a fixed value implementation is given below.
750  * @code
751  * <activityCoefficients model="Pitzer">
752  * <!-- A_Debye units = sqrt(kg/gmol) -->
753  * <A_Debye> 1.172576 </A_Debye>
754  * <!-- object description continues -->
755  * </activityCoefficients>
756  * @endcode
757  *
758  * An example of a variable value implementation within the HMWSoln object is
759  * given below. The model attribute, "water", triggers the full implementation.
760  *
761  * @code
762  * <activityCoefficients model="Pitzer">
763  * <!-- A_Debye units = sqrt(kg/gmol) -->
764  * <A_Debye model="water" />
765  * <!-- object description continues -->
766  * </activityCoefficients>
767  * @endcode
768  *
769  * ### Temperature and Pressure Dependence of the Activity Coefficients
770  *
771  * Temperature dependence of the activity coefficients leads to nonzero terms
772  * for the excess enthalpy and entropy of solution. This means that the partial
773  * molar enthalpies, entropies, and heat capacities are all non-trivial to
774  * compute. The following formulas are used.
775  *
776  * The partial molar enthalpy, \f$ \bar s_k(T,P) \f$:
777  *
778  * \f[
779  * \bar h_k(T,P) = h^{\triangle}_k(T,P)
780  * - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
781  * \f]
782  * The solvent partial molar enthalpy is equal to
783  * \f[
784  * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o)}{dT}
785  * = h^{o}_o(T,P)
786  * + R T^2 (\sum_{k \neq o} m_k) \tilde{M_o} (\frac{d \phi}{dT})
787  * \f]
788  *
789  * The partial molar entropy, \f$ \bar s_k(T,P) \f$:
790  *
791  * \f[
792  * \bar s_k(T,P) = s^{\triangle}_k(T,P)
793  * - R \ln( \gamma^{\triangle}_k \frac{m_k}{m^{\triangle}}))
794  * - R T \frac{d \ln(\gamma^{\triangle}_k) }{dT}
795  * \f]
796  * \f[
797  * \bar s_o(T,P) = s^o_o(T,P) - R \ln(a_o)
798  * - R T \frac{d \ln(a_o)}{dT}
799  * \f]
800  *
801  * The partial molar heat capacity, \f$ C_{p,k}(T,P)\f$:
802  *
803  * \f[
804  * \bar C_{p,k}(T,P) = C^{\triangle}_{p,k}(T,P)
805  * - 2 R T \frac{d \ln( \gamma^{\triangle}_k)}{dT}
806  * - R T^2 \frac{d^2 \ln(\gamma^{\triangle}_k) }{{dT}^2}
807  * \f]
808  * \f[
809  * \bar C_{p,o}(T,P) = C^o_{p,o}(T,P)
810  * - 2 R T \frac{d \ln(a_o)}{dT}
811  * - R T^2 \frac{d^2 \ln(a_o)}{{dT}^2}
812  * \f]
813  *
814  * The pressure dependence of the activity coefficients leads to non-zero terms
815  * for the excess Volume of the solution. Therefore, the partial molar volumes
816  * are functions of the pressure derivatives of the activity coefficients.
817  * \f[
818  * \bar V_k(T,P) = V^{\triangle}_k(T,P)
819  * + R T \frac{d \ln(\gamma^{\triangle}_k) }{dP}
820  * \f]
821  * \f[
822  * \bar V_o(T,P) = V^o_o(T,P)
823  * + R T \frac{d \ln(a_o)}{dP}
824  * \f]
825  *
826  * The majority of work for these functions take place in the internal routines
827  * that calculate the first and second derivatives of the log of the activity
828  * coefficients wrt temperature, s_update_dlnMolalityActCoeff_dT(),
829  * s_update_d2lnMolalityActCoeff_dT2(), and the first derivative of the log
830  * activity coefficients wrt pressure, s_update_dlnMolalityActCoeff_dP().
831  *
832  * ## %Application within Kinetics Managers
833  *
834  * For the time being, we have set the standard concentration for all solute
835  * species in this phase equal to the default concentration of the solvent at
836  * the system temperature and pressure multiplied by Mnaught (kg solvent / gmol
837  * solvent). The solvent standard concentration is just equal to its standard
838  * state concentration.
839  *
840  * This means that the kinetics operator essentially works on an generalized
841  * concentration basis (kmol / m3), with units for the kinetic rate constant
842  * specified as if all reactants (solvent or solute) are on a concentration
843  * basis (kmol /m3). The concentration will be modified by the activity
844  * coefficients.
845  *
846  * For example, a bulk-phase binary reaction between liquid solute species *j*
847  * and *k*, producing a new liquid solute species *l* would have the following
848  * equation for its rate of progress variable, \f$ R^1 \f$, which has units of
849  * kmol m-3 s-1.
850  *
851  * \f[
852  * R^1 = k^1 C_j^a C_k^a = k^1 (C^o_o \tilde{M}_o a_j) (C^o_o \tilde{M}_o a_k)
853  * \f]
854  *
855  * where
856  *
857  * \f[
858  * C_j^a = C^o_o \tilde{M}_o a_j \quad and \quad C_k^a = C^o_o \tilde{M}_o a_k
859  * \f]
860  *
861  * \f$ C_j^a \f$ is the activity concentration of species *j*, and
862  * \f$ C_k^a \f$ is the activity concentration of species *k*. \f$ C^o_o \f$ is
863  * the concentration of water at 298 K and 1 atm. \f$ \tilde{M}_o \f$ has units
864  * of kg solvent per gmol solvent and is equal to
865  *
866  * \f[
867  * \tilde{M}_o = \frac{M_o}{1000}
868  * \f]
869  *
870  * \f$ a_j \f$ is the activity of species *j* at the current temperature and
871  * pressure and concentration of the liquid phase is given by the molality based
872  * activity coefficient multiplied by the molality of the jth species.
873  *
874  * \f[
875  * a_j = \gamma_j^\triangle m_j = \gamma_j^\triangle \frac{n_j}{\tilde{M}_o n_o}
876  * \f]
877  *
878  * \f$k^1 \f$ has units of m^3/kmol/s.
879  *
880  * Therefore the generalized activity concentration of a solute species has the following form
881  *
882  * \f[
883  * C_j^a = C^o_o \frac{\gamma_j^\triangle n_j}{n_o}
884  * \f]
885  *
886  * The generalized activity concentration of the solvent has the same units, but it's a simpler form
887  *
888  * \f[
889  * C_o^a = C^o_o a_o
890  * \f]
891  *
892  * The reverse rate constant can then be obtained from the law of microscopic reversibility
893  * and the equilibrium expression for the system.
894  *
895  * \f[
896  * \frac{a_j a_k}{ a_l} = K^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} )
897  * \f]
898  *
899  * \f$ K^{o,1} \f$ is the dimensionless form of the equilibrium constant.
900  *
901  * \f[
902  * R^{-1} = k^{-1} C_l^a = k^{-1} (C_o \tilde{M}_o a_l)
903  * \f]
904  *
905  * where
906  *
907  * \f[
908  * k^{-1} = k^1 K^{o,1} C_o \tilde{M}_o
909  * \f]
910  *
911  * \f$ k^{-1} \f$ has units of 1/s.
912  *
913  * Note, this treatment may be modified in the future, as events dictate.
914  *
915  * ## Instantiation of the Class
916  *
917  * The constructor for this phase is now located in the default ThermoFactory
918  * for %Cantera. The following code snippet may be used to initialize the phase
919  * using the default construction technique within %Cantera.
920  *
921  * @code
922  * ThermoPhase *HMW = newPhase("HMW_NaCl.xml", "NaCl_electrolyte");
923  * @endcode
924  *
925  * A new HMWSoln object may be created by the following code snippets:
926  *
927  * @code
928  * HMWSoln *HMW = new HMWSoln("HMW_NaCl.xml", "NaCl_electrolyte");
929  * @endcode
930  *
931  * or
932  *
933  * @code
934  * XML_Node *xm = get_XML_NameID("phase", "HMW_NaCl.xml#NaCl_electrolyte", 0);
935  * HMWSoln *dh = new HMWSoln(*xm);
936  * @endcode
937  *
938  * or by the following call to importPhase():
939  *
940  * @code
941  * XML_Node *xm = get_XML_NameID("phase", "HMW_NaCl.xml#NaCl_electrolyte", 0);
942  * HMWSoln dhphase;
943  * importPhase(*xm, &dhphase);
944  * @endcode
945  *
946  * ## XML Example
947  *
948  * The phase model name for this is called StoichSubstance. It must be supplied
949  * as the model attribute of the thermo XML element entry. Within the phase XML
950  * block, the density of the phase must be specified. An example of an XML file
951  * this phase is given below.
952  *
953  * @code
954  * <phase id="NaCl_electrolyte" dim="3">
955  * <speciesArray datasrc="#species_waterSolution">
956  * H2O(L) Na+ Cl- H+ OH-
957  * </speciesArray>
958  * <state>
959  * <temperature units="K"> 300 </temperature>
960  * <pressure units="Pa">101325.0</pressure>
961  * <soluteMolalities>
962  * Na+:3.0
963  * Cl-:3.0
964  * H+:1.0499E-8
965  * OH-:1.3765E-6
966  * </soluteMolalities>
967  * </state>
968  * <!-- thermo model identifies the inherited class
969  * from ThermoPhase that will handle the thermodynamics.
970  * -->
971  * <thermo model="HMW">
972  * <standardConc model="solvent_volume" />
973  * <activityCoefficients model="Pitzer" TempModel="complex1">
974  * <!-- Pitzer Coefficients
975  * These coefficients are from Pitzer's main
976  * paper, in his book.
977  * -->
978  * <A_Debye model="water" />
979  * <ionicRadius default="3.042843" units="Angstroms">
980  * </ionicRadius>
981  * <binarySaltParameters cation="Na+" anion="Cl-">
982  * <beta0> 0.0765, 0.008946, -3.3158E-6,
983  * -777.03, -4.4706
984  * </beta0>
985  * <beta1> 0.2664, 6.1608E-5, 1.0715E-6 </beta1>
986  * <beta2> 0.0 </beta2>
987  * <Cphi> 0.00127, -4.655E-5, 0.0,
988  * 33.317, 0.09421
989  * </Cphi>
990  * <Alpha1> 2.0 </Alpha1>
991  * </binarySaltParameters>
992  *
993  * <binarySaltParameters cation="H+" anion="Cl-">
994  * <beta0> 0.1775, 0.0, 0.0, 0.0, 0.0</beta0>
995  * <beta1> 0.2945, 0.0, 0.0 </beta1>
996  * <beta2> 0.0 </beta2>
997  * <Cphi> 0.0008, 0.0, 0.0, 0.0, 0.0 </Cphi>
998  * <Alpha1> 2.0 </Alpha1>
999  * </binarySaltParameters>
1000  *
1001  * <binarySaltParameters cation="Na+" anion="OH-">
1002  * <beta0> 0.0864, 0.0, 0.0, 0.0, 0.0 </beta0>
1003  * <beta1> 0.253, 0.0, 0.0 </beta1>
1004  * <beta2> 0.0 </beta2>
1005  * <Cphi> 0.0044, 0.0, 0.0, 0.0, 0.0 </Cphi>
1006  * <Alpha1> 2.0 </Alpha1>
1007  * </binarySaltParameters>
1008  *
1009  * <thetaAnion anion1="Cl-" anion2="OH-">
1010  * <Theta> -0.05 </Theta>
1011  * </thetaAnion>
1012  *
1013  * <psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
1014  * <Theta> -0.05 </Theta>
1015  * <Psi> -0.006 </Psi>
1016  * </psiCommonCation>
1017  *
1018  * <thetaCation cation1="Na+" cation2="H+">
1019  * <Theta> 0.036 </Theta>
1020  * </thetaCation>
1021  *
1022  * <psiCommonAnion anion="Cl-" cation1="Na+" cation2="H+">
1023  * <Theta> 0.036 </Theta>
1024  * <Psi> -0.004 </Psi>
1025  * </psiCommonAnion>
1026  *
1027  * </activityCoefficients>
1028  *
1029  * <solvent> H2O(L) </solvent>
1030  * </thermo>
1031  * <elementArray datasrc="elements.xml"> O H Na Cl </elementArray>
1032  * <kinetics model="none" >
1033  * </kinetics>
1034  * </phase>
1035  * @endcode
1036  * @ingroup thermoprops
1037  */
1038 class HMWSoln : public MolalityVPSSTP
1039 {
1040 public:
1041  //! Default Constructor
1042  HMWSoln();
1043  ~HMWSoln();
1044 
1045  //! Construct and initialize an HMWSoln ThermoPhase object
1046  //! directly from an ASCII input file
1047  /*!
1048  * This constructor is a shell that calls the routine initThermo(), with
1049  * a reference to the XML database to get the info for the phase.
1050  *
1051  * @param inputFile Name of the input file containing the phase XML data
1052  * to set up the object
1053  * @param id ID of the phase in the input file. Defaults to the
1054  * empty string.
1055  */
1056  HMWSoln(const std::string& inputFile, const std::string& id = "");
1057 
1058  //! Construct and initialize an HMWSoln ThermoPhase object
1059  //! directly from an XML database
1060  /*!
1061  * @param phaseRef XML phase node containing the description of the phase
1062  * @param id id attribute containing the name of the phase.
1063  * (default is the empty string)
1064  */
1065  HMWSoln(XML_Node& phaseRef, const std::string& id = "");
1066 
1067  //! @name Utilities
1068  //! @{
1069 
1070  virtual std::string type() const {
1071  return "HMWSoln";
1072  }
1073 
1074  //! @}
1075  //! @name Molar Thermodynamic Properties of the Solution
1076  //! @{
1077 
1078  /// Molar enthalpy. Units: J/kmol.
1079  /**
1080  * Molar enthalpy of the solution. Units: J/kmol.
1081  * (HKM -> Bump up to Parent object)
1082  */
1083  virtual doublereal enthalpy_mole() const;
1084 
1085  /**
1086  * Excess molar enthalpy of the solution from
1087  * the mixing process. Units: J/ kmol.
1088  *
1089  * Note this is kmol of the total solution.
1090  */
1091  virtual doublereal relative_enthalpy() const;
1092 
1093  /**
1094  * Excess molar enthalpy of the solution from
1095  * the mixing process on a molality basis.
1096  * Units: J/ (kmol add salt).
1097  *
1098  * Note this is kmol of the guessed at salt composition
1099  */
1100  virtual doublereal relative_molal_enthalpy() const;
1101 
1102  /// Molar entropy. Units: J/kmol/K.
1103  /**
1104  * Molar entropy of the solution. Units: J/kmol/K. For an ideal, constant
1105  * partial molar volume solution mixture with pure species phases which
1106  * exhibit zero volume expansivity:
1107  * \f[
1108  * \hat s(T, P, X_k) = \sum_k X_k \hat s^0_k(T)
1109  * - \hat R \sum_k X_k log(X_k)
1110  * \f]
1111  * The reference-state pure-species entropies \f$ \hat s^0_k(T,p_{ref}) \f$
1112  * are computed by the species thermodynamic property manager. The pure
1113  * species entropies are independent of temperature since the volume
1114  * expansivities are equal to zero.
1115  * @see MultiSpeciesThermo
1116  *
1117  * (HKM -> Bump up to Parent object)
1118  */
1119  virtual doublereal entropy_mole() const;
1120 
1121  /// Molar Gibbs function. Units: J/kmol.
1122  /*!
1123  * (HKM -> Bump up to Parent object)
1124  */
1125  virtual doublereal gibbs_mole() const;
1126 
1127  virtual doublereal cp_mole() const;
1128 
1129  /// Molar heat capacity at constant volume. Units: J/kmol/K.
1130  /*!
1131  * (HKM -> Bump up to Parent object)
1132  */
1133  virtual doublereal cv_mole() const;
1134 
1135  //!@}
1136  //! @name Mechanical Equation of State Properties
1137  /*!
1138  * In this equation of state implementation, the density is a function
1139  * only of the mole fractions. Therefore, it can't be an independent
1140  * variable. Instead, the pressure is used as the independent variable.
1141  * Functions which try to set the thermodynamic state by calling
1142  * setDensity() may cause an exception to be thrown.
1143  */
1144  //!@{
1145 
1146 protected:
1147  /**
1148  * Calculate the density of the mixture using the partial
1149  * molar volumes and mole fractions as input
1150  *
1151  * The formula for this is
1152  *
1153  * \f[
1154  * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}}
1155  * \f]
1156  *
1157  * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are the molecular
1158  * weights, and \f$V_k\f$ are the pure species molar volumes.
1159  *
1160  * Note, the basis behind this formula is that in an ideal solution the
1161  * partial molar volumes are equal to the pure species molar volumes. We
1162  * have additionally specified in this class that the pure species molar
1163  * volumes are independent of temperature and pressure.
1164  *
1165  * NOTE: This is a non-virtual function, which is not a member of the
1166  * ThermoPhase base class.
1167  */
1168  void calcDensity();
1169 
1170 public:
1171  //! Set the internally stored density (kg/m^3) of the phase.
1172  /*!
1173  * Overridden setDensity() function is necessary because the density is not
1174  * an independent variable.
1175  *
1176  * This function will now throw an error condition.
1177  *
1178  * Note, in general, setting the phase density is now a nonlinear
1179  * calculation. P and T are the fundamental variables. This routine should
1180  * be revamped to do the nonlinear problem.
1181  *
1182  * @todo May have to adjust the strategy here to make the eos for these
1183  * materials slightly compressible, in order to create a condition where
1184  * the density is a function of the pressure.
1185  * @todo Now have a compressible ss equation for liquid water. Therefore,
1186  * this phase is compressible. May still want to change the
1187  * independent variable however.
1188  *
1189  * @param rho Input density (kg/m^3).
1190  */
1191  virtual void setDensity(const doublereal rho);
1192 
1193  //! Set the internally stored molar density (kmol/m^3) for the phase.
1194  /**
1195  * Overridden setMolarDensity() function is necessary because of the
1196  * underlying water model.
1197  *
1198  * This function will now throw an error condition if the input isn't
1199  * exactly equal to the current molar density.
1200  *
1201  * @param conc Input molar density (kmol/m^3).
1202  */
1203  virtual void setMolarDensity(const doublereal conc);
1204 
1205  /**
1206  * @}
1207  * @name Activities, Standard States, and Activity Concentrations
1208  *
1209  * The activity \f$a_k\f$ of a species in solution is related to the
1210  * chemical potential by \f[ \mu_k = \mu_k^0(T) + \hat R T \log a_k. \f] The
1211  * quantity \f$\mu_k^0(T,P)\f$ is the chemical potential at unit activity,
1212  * which depends only on temperature and the pressure. Activity is assumed
1213  * to be molality-based here.
1214  * @{
1215  */
1216 
1217  //! This method returns an array of generalized activity concentrations
1218  /*!
1219  * The generalized activity concentrations, \f$ C_k^a\f$, are defined such
1220  * that \f$ a_k = C^a_k / C^0_k, \f$ where \f$ C^0_k \f$ is a standard
1221  * concentration defined below. These generalized concentrations are used
1222  * by kinetics manager classes to compute the forward and reverse rates of
1223  * elementary reactions.
1224  *
1225  * The generalized activity concentration of a solute species has the
1226  * following form
1227  *
1228  * \f[
1229  * C_j^a = C^o_o \frac{\gamma_j^\triangle n_j}{n_o}
1230  * \f]
1231  *
1232  * The generalized activity concentration of the solvent has the same units,
1233  * but it's a simpler form
1234  *
1235  * \f[
1236  * C_o^a = C^o_o a_o
1237  * \f]
1238  *
1239  * @param c Array of generalized concentrations. The
1240  * units are kmol m-3 for both the solvent and the solute species
1241  */
1242  virtual void getActivityConcentrations(doublereal* c) const;
1243 
1244  //! Return the standard concentration for the kth species
1245  /*!
1246  * The standard concentration \f$ C^0_k \f$ used to normalize the activity
1247  * (i.e., generalized) concentration for use
1248  *
1249  * We have set the standard concentration for all solute species in this
1250  * phase equal to the default concentration of the solvent at the system
1251  * temperature and pressure multiplied by Mnaught (kg solvent / gmol
1252  * solvent). The solvent standard concentration is just equal to its
1253  * standard state concentration.
1254  *
1255  * \f[
1256  * C_j^0 = C^o_o \tilde{M}_o \quad and C_o^0 = C^o_o
1257  * \f]
1258  *
1259  * The consequence of this is that the standard concentrations have unequal
1260  * units between the solvent and the solute. However, both the solvent and
1261  * the solute activity concentrations will have the same units of kmol/kg^3.
1262  *
1263  * This means that the kinetics operator essentially works on an generalized
1264  * concentration basis (kmol / m3), with units for the kinetic rate constant
1265  * specified as if all reactants (solvent or solute) are on a concentration
1266  * basis (kmol /m3). The concentration will be modified by the activity
1267  * coefficients.
1268  *
1269  * For example, a bulk-phase binary reaction between liquid solute species
1270  * *j* and *k*, producing a new liquid solute species *l* would have the
1271  * following equation for its rate of progress variable, \f$ R^1 \f$, which
1272  * has units of kmol m-3 s-1.
1273  *
1274  * \f[
1275  * R^1 = k^1 C_j^a C_k^a = k^1 (C^o_o \tilde{M}_o a_j) (C^o_o \tilde{M}_o a_k)
1276  * \f]
1277  *
1278  * where
1279  *
1280  * \f[
1281  * C_j^a = C^o_o \tilde{M}_o a_j \quad and \quad C_k^a = C^o_o \tilde{M}_o a_k
1282  * \f]
1283  *
1284  * \f$ C_j^a \f$ is the activity concentration of species *j*, and
1285  * \f$ C_k^a \f$ is the activity concentration of species *k*. \f$ C^o_o \f$
1286  * is the concentration of water at 298 K and 1 atm. \f$ \tilde{M}_o \f$ has
1287  * units of kg solvent per gmol solvent and is equal to
1288  *
1289  * \f[
1290  * \tilde{M}_o = \frac{M_o}{1000}
1291  * \f]
1292  *
1293  * \f$ a_j \f$ is
1294  * the activity of species *j* at the current temperature and pressure
1295  * and concentration of the liquid phase is given by the molality based
1296  * activity coefficient multiplied by the molality of the jth species.
1297  *
1298  * \f[
1299  * a_j = \gamma_j^\triangle m_j = \gamma_j^\triangle \frac{n_j}{\tilde{M}_o n_o}
1300  * \f]
1301  *
1302  * \f$k^1 \f$ has units of m^3/kmol/s.
1303  *
1304  * Therefore the generalized activity concentration of a solute species has
1305  * the following form
1306  *
1307  * \f[
1308  * C_j^a = C^o_o \frac{\gamma_j^\triangle n_j}{n_o}
1309  * \f]
1310  *
1311  * The generalized activity concentration of the solvent has the same units,
1312  * but it's a simpler form
1313  *
1314  * \f[
1315  * C_o^a = C^o_o a_o
1316  * \f]
1317  *
1318  * @param k Optional parameter indicating the species. The default is to
1319  * assume this refers to species 0.
1320  * @returns the standard Concentration in units of m^3/kmol.
1321  *
1322  * @param k Species index
1323  */
1324  virtual doublereal standardConcentration(size_t k=0) const;
1325 
1326  //! Get the array of non-dimensional activities at the current solution
1327  //! temperature, pressure, and solution concentration.
1328  /*!
1329  *
1330  * We resolve this function at this level by calling on the
1331  * activityConcentration function. However, derived classes may want to
1332  * override this default implementation.
1333  *
1334  * (note solvent is on molar scale).
1335  *
1336  * @param ac Output vector of activities. Length: m_kk.
1337  */
1338  virtual void getActivities(doublereal* ac) const;
1339 
1340  //! @}
1341  //! @name Partial Molar Properties of the Solution
1342  //! @{
1343 
1344  //! Get the species chemical potentials. Units: J/kmol.
1345  /*!
1346  *
1347  * This function returns a vector of chemical potentials of the
1348  * species in solution.
1349  *
1350  * \f[
1351  * \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} m_k)
1352  * \f]
1353  *
1354  * @param mu Output vector of species chemical
1355  * potentials. Length: m_kk. Units: J/kmol
1356  */
1357  virtual void getChemPotentials(doublereal* mu) const;
1358 
1359  //! Returns an array of partial molar enthalpies for the species
1360  //! in the mixture. Units (J/kmol)
1361  /*!
1362  * For this phase, the partial molar enthalpies are equal to the standard
1363  * state enthalpies modified by the derivative of the molality-based
1364  * activity coefficient wrt temperature
1365  *
1366  * \f[
1367  * \bar h_k(T,P) = h^{\triangle}_k(T,P)
1368  * - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
1369  * \f]
1370  * The solvent partial molar enthalpy is equal to
1371  * \f[
1372  * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o)}{dT}
1373  * = h^{o}_o(T,P)
1374  * + R T^2 (\sum_{k \neq o} m_k) \tilde{M_o} (\frac{d \phi}{dT})
1375  * \f]
1376  *
1377  * @param hbar Output vector of species partial molar enthalpies.
1378  * Length: m_kk. units are J/kmol.
1379  */
1380  virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
1381 
1382  //! Returns an array of partial molar entropies of the species in the
1383  //! solution. Units: J/kmol/K.
1384  /*!
1385  * Maxwell's equations provide an answer for how calculate this
1386  * (p.215 Smith and Van Ness)
1387  *
1388  * d(chemPot_i)/dT = -sbar_i
1389  *
1390  * For this phase, the partial molar entropies are equal to the SS species
1391  * entropies plus the ideal solution contribution plus complicated functions
1392  * of the temperature derivative of the activity coefficients.
1393  *
1394  * \f[
1395  * \bar s_k(T,P) = s^{\triangle}_k(T,P)
1396  * - R \ln( \gamma^{\triangle}_k \frac{m_k}{m^{\triangle}}))
1397  * - R T \frac{d \ln(\gamma^{\triangle}_k) }{dT}
1398  * \f]
1399  * \f[
1400  * \bar s_o(T,P) = s^o_o(T,P) - R \ln(a_o)
1401  * - R T \frac{d \ln(a_o)}{dT}
1402  * \f]
1403  *
1404  * @param sbar Output vector of species partial molar entropies.
1405  * Length = m_kk. units are J/kmol/K.
1406  */
1407  virtual void getPartialMolarEntropies(doublereal* sbar) const;
1408 
1409  //! Return an array of partial molar volumes for the species in the mixture.
1410  //! Units: m^3/kmol.
1411  /*!
1412  * For this solution, the partial molar volumes are functions of the
1413  * pressure derivatives of the activity coefficients.
1414  *
1415  * \f[
1416  * \bar V_k(T,P) = V^{\triangle}_k(T,P)
1417  * + R T \frac{d \ln(\gamma^{\triangle}_k) }{dP}
1418  * \f]
1419  * \f[
1420  * \bar V_o(T,P) = V^o_o(T,P)
1421  * + R T \frac{d \ln(a_o)}{dP}
1422  * \f]
1423  *
1424  * @param vbar Output vector of species partial molar volumes.
1425  * Length = m_kk. units are m^3/kmol.
1426  */
1427  virtual void getPartialMolarVolumes(doublereal* vbar) const;
1428 
1429  //! Return an array of partial molar heat capacities for the species in the
1430  //! mixture. Units: J/kmol/K
1431  /*!
1432  * The following formulas are implemented within the code.
1433  *
1434  * \f[
1435  * \bar C_{p,k}(T,P) = C^{\triangle}_{p,k}(T,P)
1436  * - 2 R T \frac{d \ln( \gamma^{\triangle}_k)}{dT}
1437  * - R T^2 \frac{d^2 \ln(\gamma^{\triangle}_k) }{{dT}^2}
1438  * \f]
1439  * \f[
1440  * \bar C_{p,o}(T,P) = C^o_{p,o}(T,P)
1441  * - 2 R T \frac{d \ln(a_o)}{dT}
1442  * - R T^2 \frac{d^2 \ln(a_o)}{{dT}^2}
1443  * \f]
1444  *
1445  * @param cpbar Output vector of species partial molar heat capacities at
1446  * constant pressure. Length = m_kk. units are J/kmol/K.
1447  */
1448  virtual void getPartialMolarCp(doublereal* cpbar) const;
1449 
1450 public:
1451  //@}
1452 
1453  //! Get the saturation pressure for a given temperature.
1454  /*!
1455  * Note the limitations of this function. Stability considerations
1456  * concerning multiphase equilibrium are ignored in this calculation.
1457  * Therefore, the call is made directly to the SS of water underneath. The
1458  * object is put back into its original state at the end of the call.
1459  *
1460  * @todo This is probably not implemented correctly. The stability of the
1461  * salt should be added into this calculation. The underlying water
1462  * model may be called to get the stability of the pure water
1463  * solution, if needed.
1464  *
1465  * @param T Temperature (kelvin)
1466  */
1467  virtual doublereal satPressure(doublereal T);
1468 
1469  /*
1470  * -------------- Utilities -------------------------------
1471  */
1472 
1473  void setBinarySalt(const std::string& sp1, const std::string& sp2,
1474  size_t nParams, double* beta0, double* beta1, double* beta2,
1475  double* Cphi, double alpha1, double alpha2);
1476  void setTheta(const std::string& sp1, const std::string& sp2,
1477  size_t nParams, double* theta);
1478  void setPsi(const std::string& sp1, const std::string& sp2,
1479  const std::string& sp3, size_t nParams, double* psi);
1480  void setLambda(const std::string& sp1, const std::string& sp2,
1481  size_t nParams, double* lambda);
1482  void setMunnn(const std::string& sp, size_t nParams, double* munnn);
1483  void setZeta(const std::string& sp1, const std::string& sp2,
1484  const std::string& sp3, size_t nParams, double* psi);
1485 
1486  void setPitzerTempModel(const std::string& model);
1487  void setPitzerRefTemperature(double Tref) {
1488  m_TempPitzerRef = Tref;
1489  }
1490 
1491  //! Set the A_Debye parameter. If a negative value is provided, enables
1492  //! calculation of A_Debye using the detailed water equation of state.
1493  void setA_Debye(double A);
1494 
1495  void setMaxIonicStrength(double Imax) {
1496  m_maxIionicStrength = Imax;
1497  }
1498 
1499  void setCroppingCoefficients(double ln_gamma_k_min, double ln_gamma_k_max,
1500  double ln_gamma_o_min, double ln_gamma_o_max);
1501 
1502  virtual void initThermo();
1503 
1504  //! Initialize the phase parameters from an XML file.
1505  /*!
1506  * This gets called from importPhase(). It processes the XML file after the
1507  * species are set up. This is the main routine for reading in activity
1508  * coefficient parameters.
1509  *
1510  * @param phaseNode This object must be the phase node of a complete XML
1511  * tree description of the phase, including all of the species
1512  * data. In other words while "phase" must point to an XML phase
1513  * object, it must have sibling nodes "speciesData" that
1514  * describe the species in the phase.
1515  * @param id ID of the phase. If nonnull, a check is done to see if
1516  * phaseNode is pointing to the phase with the correct id.
1517  */
1518  virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
1519 
1520  //! Value of the Debye Huckel constant as a function of temperature
1521  //! and pressure.
1522  /*!
1523  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1524  *
1525  * Units = sqrt(kg/gmol)
1526  *
1527  * @param temperature Temperature of the derivative calculation
1528  * or -1 to indicate the current temperature
1529  * @param pressure Pressure of the derivative calculation
1530  * or -1 to indicate the current pressure
1531  */
1532  virtual double A_Debye_TP(double temperature = -1.0,
1533  double pressure = -1.0) const;
1534 
1535  //! Value of the derivative of the Debye Huckel constant with respect to
1536  //! temperature as a function of temperature and pressure.
1537  /*!
1538  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1539  *
1540  * Units = sqrt(kg/gmol)
1541  *
1542  * @param temperature Temperature of the derivative calculation
1543  * or -1 to indicate the current temperature
1544  * @param pressure Pressure of the derivative calculation
1545  * or -1 to indicate the current pressure
1546  */
1547  virtual double dA_DebyedT_TP(double temperature = -1.0,
1548  double pressure = -1.0) const;
1549 
1550  /**
1551  * Value of the derivative of the Debye Huckel constant with respect to
1552  * pressure, as a function of temperature and pressure.
1553  *
1554  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1555  *
1556  * Units = sqrt(kg/gmol)
1557  *
1558  * @param temperature Temperature of the derivative calculation
1559  * or -1 to indicate the current temperature
1560  * @param pressure Pressure of the derivative calculation
1561  * or -1 to indicate the current pressure
1562  */
1563  virtual double dA_DebyedP_TP(double temperature = -1.0,
1564  double pressure = -1.0) const;
1565 
1566  /**
1567  * Return Pitzer's definition of A_L. This is basically the
1568  * derivative of the A_phi multiplied by 4 R T**2
1569  *
1570  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1571  * dA_phidT = d(A_Debye)/dT / 3.0
1572  * A_L = dA_phidT * (4 * R * T * T)
1573  *
1574  * Units = sqrt(kg/gmol) (RT)
1575  *
1576  * @param temperature Temperature of the derivative calculation
1577  * or -1 to indicate the current temperature
1578  * @param pressure Pressure of the derivative calculation
1579  * or -1 to indicate the current pressure
1580  */
1581  double ADebye_L(double temperature = -1.0,
1582  double pressure = -1.0) const;
1583 
1584  /**
1585  * Return Pitzer's definition of A_J. This is basically the temperature
1586  * derivative of A_L, and the second derivative of A_phi
1587  *
1588  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1589  * dA_phidT = d(A_Debye)/dT / 3.0
1590  * A_J = 2 A_L/T + 4 * R * T * T * d2(A_phi)/dT2
1591  *
1592  * Units = sqrt(kg/gmol) (R)
1593  *
1594  * @param temperature Temperature of the derivative calculation
1595  * or -1 to indicate the current temperature
1596  * @param pressure Pressure of the derivative calculation
1597  * or -1 to indicate the current pressure
1598  */
1599  double ADebye_J(double temperature = -1.0,
1600  double pressure = -1.0) const;
1601 
1602  /**
1603  * Return Pitzer's definition of A_V. This is the derivative wrt pressure of
1604  * A_phi multiplied by - 4 R T
1605  *
1606  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1607  * dA_phidT = d(A_Debye)/dP / 3.0
1608  * A_V = - dA_phidP * (4 * R * T)
1609  *
1610  * Units = sqrt(kg/gmol) (RT) / Pascal
1611  *
1612  * @param temperature Temperature of the derivative calculation
1613  * or -1 to indicate the current temperature
1614  * @param pressure Pressure of the derivative calculation
1615  * or -1 to indicate the current pressure
1616  */
1617  double ADebye_V(double temperature = -1.0,
1618  double pressure = -1.0) const;
1619 
1620  //! Value of the 2nd derivative of the Debye Huckel constant with respect to
1621  //! temperature as a function of temperature and pressure.
1622  /*!
1623  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1624  *
1625  * Units = sqrt(kg/gmol)
1626  *
1627  * @param temperature Temperature of the derivative calculation
1628  * or -1 to indicate the current temperature
1629  * @param pressure Pressure of the derivative calculation
1630  * or -1 to indicate the current pressure
1631  */
1632  virtual double d2A_DebyedT2_TP(double temperature = -1.0,
1633  double pressure = -1.0) const;
1634 
1635  //! Print out all of the input Pitzer coefficients.
1636  void printCoeffs() const;
1637 
1638  //! Get the array of unscaled non-dimensional molality based
1639  //! activity coefficients at the current solution temperature,
1640  //! pressure, and solution concentration.
1641  /*!
1642  * See Denbigh p. 278 for a thorough discussion. This class must be
1643  * overridden in classes which derive from MolalityVPSSTP. This function
1644  * takes over from the molar-based activity coefficient calculation,
1645  * getActivityCoefficients(), in derived classes.
1646  *
1647  * @param acMolality Output vector containing the molality based activity coefficients.
1648  * length: m_kk.
1649  */
1650  void getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const;
1651 
1652 private:
1653  //! Apply the current phScale to a set of activity Coefficients
1654  /*!
1655  * See the Eq3/6 Manual for a thorough discussion.
1656  */
1657  void s_updateScaling_pHScaling() const;
1658 
1659  //! Apply the current phScale to a set of derivatives of the activity
1660  //! Coefficients wrt temperature
1661  /*!
1662  * See the Eq3/6 Manual for a thorough discussion of the need
1663  */
1664  void s_updateScaling_pHScaling_dT() const;
1665 
1666  //! Apply the current phScale to a set of 2nd derivatives of the activity
1667  //! Coefficients wrt temperature
1668  /*!
1669  * See the Eq3/6 Manual for a thorough discussion of the need
1670  */
1671  void s_updateScaling_pHScaling_dT2() const;
1672 
1673  //! Apply the current phScale to a set of derivatives of the activity
1674  //! Coefficients wrt pressure
1675  /*!
1676  * See the Eq3/6 Manual for a thorough discussion of the need
1677  */
1678  void s_updateScaling_pHScaling_dP() const;
1679 
1680  //! Calculate the Chlorine activity coefficient on the NBS scale
1681  /*!
1682  * We assume here that the m_IionicMolality variable is up to date.
1683  */
1684  doublereal s_NBS_CLM_lnMolalityActCoeff() const;
1685 
1686  //! Calculate the temperature derivative of the Chlorine activity
1687  //! coefficient on the NBS scale
1688  /*!
1689  * We assume here that the m_IionicMolality variable is up to date.
1690  */
1691  doublereal s_NBS_CLM_dlnMolalityActCoeff_dT() const;
1692 
1693  //! Calculate the second temperature derivative of the Chlorine activity
1694  //! coefficient on the NBS scale
1695  /*!
1696  * We assume here that the m_IionicMolality variable is up to date.
1697  */
1698  doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const;
1699 
1700  //! Calculate the pressure derivative of the Chlorine activity coefficient
1701  /*!
1702  * We assume here that the m_IionicMolality variable is up to date.
1703  */
1704  doublereal s_NBS_CLM_dlnMolalityActCoeff_dP() const;
1705 
1706  //@}
1707 
1708 private:
1709  /**
1710  * This is the form of the temperature dependence of Pitzer parameterization
1711  * used in the model.
1712  *
1713  * PITZER_TEMP_CONSTANT 0
1714  * PITZER_TEMP_LINEAR 1
1715  * PITZER_TEMP_COMPLEX1 2
1716  */
1718 
1719  //! Current value of the ionic strength on the molality scale Associated
1720  //! Salts, if present in the mechanism, don't contribute to the value of the
1721  //! ionic strength in this version of the Ionic strength.
1722  mutable double m_IionicMolality;
1723 
1724  //! Maximum value of the ionic strength allowed in the calculation of the
1725  //! activity coefficients.
1727 
1728  //! Reference Temperature for the Pitzer formulations.
1730 
1731 public:
1732  /**
1733  * Form of the constant outside the Debye-Huckel term called A. It's
1734  * normally a function of temperature and pressure. However, it can be set
1735  * from the input file in order to aid in numerical comparisons. Acceptable
1736  * forms:
1737  *
1738  * A_DEBYE_CONST 0
1739  * A_DEBYE_WATER 1
1740  *
1741  * The A_DEBYE_WATER form may be used for water solvents with needs to cover
1742  * varying temperatures and pressures. Note, the dielectric constant of
1743  * water is a relatively strong function of T, and its variability must be
1744  * accounted for,
1745  */
1747 
1748 private:
1749  /**
1750  * A_Debye: this expression appears on the top of the ln actCoeff term in
1751  * the general Debye-Huckel expression It depends on temperature.
1752  * And, therefore, most be recalculated whenever T or P changes.
1753  * This variable is a local copy of the calculation.
1754  *
1755  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1756  *
1757  * where B_Debye = F / sqrt(epsilon R T/2)
1758  * (dw/1000)^(1/2)
1759  *
1760  * A_Debye = (1/ (8 Pi)) (2 Na * dw/1000)^(1/2)
1761  * (e * e / (epsilon * kb * T))^(3/2)
1762  *
1763  * Units = sqrt(kg/gmol)
1764  *
1765  * Nominal value = 1.172576 sqrt(kg/gmol)
1766  * based on:
1767  * epsilon/epsilon_0 = 78.54
1768  * (water at 25C)
1769  * epsilon_0 = 8.854187817E-12 C2 N-1 m-2
1770  * e = 1.60217653 E-19 C
1771  * F = 9.6485309E7 C kmol-1
1772  * R = 8.314472E3 kg m2 s-2 kmol-1 K-1
1773  * T = 298.15 K
1774  * B_Debye = 3.28640E9 sqrt(kg/gmol)/m
1775  * dw = C_0 * M_0 (density of water) (kg/m3)
1776  * = 1.0E3 at 25C
1777  */
1778  mutable double m_A_Debye;
1779 
1780  //! Water standard state calculator
1781  /*!
1782  * derived from the equation of state for water.
1783  */
1785 
1786  //! Pointer to the water property calculator
1787  std::unique_ptr<WaterProps> m_waterProps;
1788 
1789  //! vector of size m_kk, used as a temporary holding area.
1791 
1792  /**
1793  * Array of 2D data used in the Pitzer/HMW formulation. Beta0_ij[i][j] is
1794  * the value of the Beta0 coefficient for the ij salt. It will be nonzero
1795  * iff i and j are both charged and have opposite sign. The array is also
1796  * symmetric. counterIJ where counterIJ = m_counterIJ[i][j] is used to
1797  * access this array.
1798  */
1800 
1801  //! Derivative of Beta0_ij[i][j] wrt T. Vector index is counterIJ
1803 
1804  //! Derivative of Beta0_ij[i][j] wrt TT. Vector index is counterIJ
1806 
1807  //! Derivative of Beta0_ij[i][j] wrt P. Vector index is counterIJ
1809 
1810  //! Array of coefficients for Beta0, a variable in Pitzer's papers
1811  /*!
1812  * Column index is counterIJ. m_Beta0MX_ij_coeff.ptrColumn(counterIJ) is a
1813  * double* containing the vector of coefficients for the counterIJ
1814  * interaction.
1815  */
1817 
1818  //! Array of 2D data used in the Pitzer/HMW formulation. Beta1_ij[i][j] is
1819  //! the value of the Beta1 coefficient for the ij salt. It will be nonzero
1820  //! iff i and j are both charged and have opposite sign. The array is also
1821  //! symmetric. counterIJ where counterIJ = m_counterIJ[i][j] is used to
1822  //! access this array.
1824 
1825  //! Derivative of Beta1_ij[i][j] wrt T. Vector index is counterIJ
1827 
1828  //! Derivative of Beta1_ij[i][j] wrt TT. Vector index is counterIJ
1830 
1831  //! Derivative of Beta1_ij[i][j] wrt P. Vector index is counterIJ
1833 
1834  //! Array of coefficients for Beta1, a variable in Pitzer's papers
1835  /*!
1836  * Column index is counterIJ. m_Beta1MX_ij_coeff.ptrColumn(counterIJ) is a
1837  * double* containing the vector of coefficients for the counterIJ
1838  * interaction.
1839  */
1841 
1842  //! Array of 2D data used in the Pitzer/HMW formulation. Beta2_ij[i][j] is
1843  //! the value of the Beta2 coefficient for the ij salt. It will be nonzero
1844  //! iff i and j are both charged and have opposite sign, and i and j both
1845  //! have charges of 2 or more. The array is also symmetric. counterIJ where
1846  //! counterIJ = m_counterIJ[i][j] is used to access this array.
1848 
1849  //! Derivative of Beta2_ij[i][j] wrt T. Vector index is counterIJ
1851 
1852  //! Derivative of Beta2_ij[i][j] wrt TT. Vector index is counterIJ
1854 
1855  //! Derivative of Beta2_ij[i][j] wrt P. Vector index is counterIJ
1857 
1858  //! Array of coefficients for Beta2, a variable in Pitzer's papers
1859  /*!
1860  * column index is counterIJ. m_Beta2MX_ij_coeff.ptrColumn(counterIJ) is a
1861  * double* containing the vector of coefficients for the counterIJ
1862  * interaction. This was added for the YMP database version of the code
1863  * since it contains temperature-dependent parameters for some 2-2
1864  * electrolytes.
1865  */
1867 
1868  // Array of 2D data used in the Pitzer/HMW formulation. Alpha1MX_ij[i][j] is
1869  // the value of the alpha1 coefficient for the ij interaction. It will be
1870  // nonzero iff i and j are both charged and have opposite sign. It is
1871  // symmetric wrt i, j. counterIJ where counterIJ = m_counterIJ[i][j] is used
1872  // to access this array.
1873  vector_fp m_Alpha1MX_ij;
1874 
1875  //! Array of 2D data used in the Pitzer/HMW formulation. Alpha2MX_ij[i][j]
1876  //! is the value of the alpha2 coefficient for the ij interaction. It will
1877  //! be nonzero iff i and j are both charged and have opposite sign, and i
1878  //! and j both have charges of 2 or more, usually. It is symmetric wrt i, j.
1879  //! counterIJ, where counterIJ = m_counterIJ[i][j], is used to access this
1880  //! array.
1882 
1883  //! Array of 2D data used in the Pitzer/HMW formulation. CphiMX_ij[i][j] is
1884  //! the value of the Cphi coefficient for the ij interaction. It will be
1885  //! nonzero iff i and j are both charged and have opposite sign, and i and j
1886  //! both have charges of 2 or more. The array is also symmetric. counterIJ
1887  //! where counterIJ = m_counterIJ[i][j] is used to access this array.
1889 
1890  //! Derivative of Cphi_ij[i][j] wrt T. Vector index is counterIJ
1892 
1893  //! Derivative of Cphi_ij[i][j] wrt TT. Vector index is counterIJ
1895 
1896  //! Derivative of Cphi_ij[i][j] wrt P. Vector index is counterIJ
1898 
1899  //! Array of coefficients for CphiMX, a parameter in the activity
1900  //! coefficient formulation
1901  /*!
1902  * Column index is counterIJ. m_CphiMX_ij_coeff.ptrColumn(counterIJ) is a
1903  * double* containing the vector of coefficients for the counterIJ
1904  * interaction.
1905  */
1907 
1908  //! Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
1909  /*!
1910  * Array of 2D data used in the Pitzer/HMW formulation. Theta_ij[i][j] is
1911  * the value of the theta coefficient for the ij interaction. It will be
1912  * nonzero for charged ions with the same sign. It is symmetric. counterIJ
1913  * where counterIJ = m_counterIJ[i][j] is used to access this array.
1914  *
1915  * HKM Recent Pitzer papers have used a functional form for Theta_ij, which
1916  * depends on the ionic strength.
1917  */
1919 
1920  //! Derivative of Theta_ij[i][j] wrt T. Vector index is counterIJ
1922 
1923  //! Derivative of Theta_ij[i][j] wrt TT. Vector index is counterIJ
1925 
1926  //! Derivative of Theta_ij[i][j] wrt P. Vector index is counterIJ
1928 
1929  //! Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
1930  /*!
1931  * Theta_ij[i][j] is the value of the theta coefficient for the ij
1932  * interaction. It will be nonzero for charged ions with the same sign. It
1933  * is symmetric. Column index is counterIJ. counterIJ where counterIJ =
1934  * m_counterIJ[i][j] is used to access this array.
1935  *
1936  * m_Theta_ij_coeff.ptrColumn(counterIJ) is a double* containing
1937  * the vector of coefficients for the counterIJ interaction.
1938  */
1940 
1941  //! Array of 3D data used in the Pitzer/HMW formulation.
1942  /*!
1943  * Psi_ijk[n] is the value of the psi coefficient for the
1944  * ijk interaction where
1945  *
1946  * n = k + j * m_kk + i * m_kk * m_kk;
1947  *
1948  * It is potentially nonzero everywhere. The first two coordinates are
1949  * symmetric wrt cations, and the last two coordinates are symmetric wrt
1950  * anions.
1951  */
1953 
1954  //! Derivative of Psi_ijk[n] wrt T. See m_Psi_ijk for reference on the
1955  //! indexing into this variable.
1957 
1958  //! Derivative of Psi_ijk[n] wrt TT. See m_Psi_ijk for reference on the
1959  //! indexing into this variable.
1961 
1962  //! Derivative of Psi_ijk[n] wrt P. See m_Psi_ijk for reference on the
1963  //! indexing into this variable.
1965 
1966  //! Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
1967  /*!
1968  * Psi_ijk[n] is the value of the psi coefficient for the
1969  * ijk interaction where
1970  *
1971  * n = k + j * m_kk + i * m_kk * m_kk;
1972  *
1973  * It is potentially nonzero everywhere. The first two coordinates are
1974  * symmetric wrt cations, and the last two coordinates are symmetric wrt
1975  * anions.
1976  *
1977  * m_Psi_ijk_coeff.ptrColumn(n) is a double* containing the vector of
1978  * coefficients for the n interaction.
1979  */
1981 
1982  //! Lambda coefficient for the ij interaction
1983  /*!
1984  * Array of 2D data used in the Pitzer/HMW formulation. Lambda_nj[n][j]
1985  * represents the lambda coefficient for the ij interaction. This is a
1986  * general interaction representing neutral species. The neutral species
1987  * occupy the first index, i.e., n. The charged species occupy the j
1988  * coordinate. neutral, neutral interactions are also included here.
1989  */
1991 
1992  //! Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij
1994 
1995  //! Derivative of Lambda_nj[i][j] wrt TT
1997 
1998  //! Derivative of Lambda_nj[i][j] wrt P
2000 
2001  //! Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
2002  /*!
2003  * Array of 2D data used in the Pitzer/HMW formulation. Lambda_ij[i][j]
2004  * represents the lambda coefficient for the ij interaction. This is a
2005  * general interaction representing neutral species. The neutral species
2006  * occupy the first index, i.e., i. The charged species occupy the j
2007  * coordinate. Neutral, neutral interactions are also included here.
2008  *
2009  * n = j + m_kk * i
2010  *
2011  * m_Lambda_ij_coeff.ptrColumn(n) is a double* containing the vector of
2012  * coefficients for the (i,j) interaction.
2013  */
2015 
2016  //! Mu coefficient for the self-ternary neutral coefficient
2017  /*!
2018  * Array of 2D data used in the Pitzer/HMW formulation. Mu_nnn[i] represents
2019  * the Mu coefficient for the nnn interaction. This is a general interaction
2020  * representing neutral species interacting with itself.
2021  */
2023 
2024  //! Mu coefficient temperature derivative for the self-ternary neutral
2025  //! coefficient
2026  /*!
2027  * Array of 2D data used in the Pitzer/HMW formulation. Mu_nnn_L[i]
2028  * represents the Mu coefficient temperature derivative for the nnn
2029  * interaction. This is a general interaction representing neutral species
2030  * interacting with itself.
2031  */
2033 
2034  //! Mu coefficient 2nd temperature derivative for the self-ternary neutral
2035  //! coefficient
2036  /*!
2037  * Array of 2D data used in the Pitzer/HMW formulation. Mu_nnn_L[i]
2038  * represents the Mu coefficient 2nd temperature derivative for the nnn
2039  * interaction. This is a general interaction representing neutral species
2040  * interacting with itself.
2041  */
2043 
2044  //! Mu coefficient pressure derivative for the self-ternary neutral
2045  //! coefficient
2046  /*!
2047  * Array of 2D data used in the Pitzer/HMW formulation. Mu_nnn_L[i]
2048  * represents the Mu coefficient pressure derivative for the nnn
2049  * interaction. This is a general interaction representing neutral species
2050  * interacting with itself.
2051  */
2053 
2054  //! Array of coefficients form_Mu_nnn term
2056 
2057  //! Logarithm of the activity coefficients on the molality scale.
2058  /*!
2059  * mutable because we change this if the composition or temperature or
2060  * pressure changes. Index is the species index
2061  */
2063 
2064  //! Logarithm of the activity coefficients on the molality scale.
2065  /*!
2066  * mutable because we change this if the composition or temperature or
2067  * pressure changes. Index is the species index
2068  */
2070 
2071  //! Derivative of the Logarithm of the activity coefficients on the molality
2072  //! scale wrt T. Index is the species index
2074 
2075  //! Derivative of the Logarithm of the activity coefficients on the molality
2076  //! scale wrt T. Index is the species index
2078 
2079  //! Derivative of the Logarithm of the activity coefficients on the molality
2080  //! scale wrt TT. Index is the species index.
2082 
2083  //! Derivative of the Logarithm of the activity coefficients on the molality
2084  //! scale wrt TT. Index is the species index
2086 
2087  //! Derivative of the Logarithm of the activity coefficients on the
2088  //! molality scale wrt P. Index is the species index
2090 
2091  //! Derivative of the Logarithm of the activity coefficients on the
2092  //! molality scale wrt P. Index is the species index
2094 
2095  // -------- Temporary Variables Used in the Activity Coeff Calc
2096 
2097  //! Cropped and modified values of the molalities used in activity
2098  //! coefficient calculations
2100 
2101  //! Boolean indicating whether the molalities are cropped or are modified
2103 
2104  //! a counter variable for keeping track of symmetric binary
2105  //! interactions amongst the solute species.
2106  /*!
2107  * n = m_kk*i + j
2108  * m_CounterIJ[n] = counterIJ
2109  */
2111 
2112  //! This is elambda, MEC
2113  mutable double elambda[17];
2114 
2115  //! This is elambda1, MEC
2116  mutable double elambda1[17];
2117 
2118  /**
2119  * Various temporary arrays used in the calculation of the Pitzer activity
2120  * coefficients. The subscript, L, denotes the same quantity's derivative
2121  * wrt temperature
2122  */
2123 
2124  //! This is the value of g(x) in Pitzer's papers. Vector index is counterIJ
2126 
2127  //! This is the value of g2(x2) in Pitzer's papers. Vector index is counterIJ
2129 
2130  //! hfunc, was called gprime in Pitzer's paper. However, it's not the
2131  //! derivative of gfunc(x), so I renamed it. Vector index is counterIJ
2133 
2134  //! hfunc2, was called gprime in Pitzer's paper. However, it's not the
2135  //! derivative of gfunc(x), so I renamed it. Vector index is counterIJ
2137 
2138  //! Intermediate variable called BMX in Pitzer's paper. This is the basic
2139  //! cation - anion interaction. Vector index is counterIJ
2141 
2142  //! Derivative of BMX_IJ wrt T. Vector index is counterIJ
2144 
2145  //! Derivative of BMX_IJ wrt TT. Vector index is counterIJ
2147 
2148  //! Derivative of BMX_IJ wrt P. Vector index is counterIJ
2150 
2151  //! Intermediate variable called BprimeMX in Pitzer's paper. Vector index is
2152  //! counterIJ
2154 
2155  //! Derivative of BprimeMX wrt T. Vector index is counterIJ
2157 
2158  //! Derivative of BprimeMX wrt TT. Vector index is counterIJ
2160 
2161  //! Derivative of BprimeMX wrt P. Vector index is counterIJ
2163 
2164  //! Intermediate variable called BphiMX in Pitzer's paper. Vector index is
2165  //! counterIJ
2167 
2168  //! Derivative of BphiMX_IJ wrt T. Vector index is counterIJ
2170 
2171  //! Derivative of BphiMX_IJ wrt TT. Vector index is counterIJ
2173 
2174  //! Derivative of BphiMX_IJ wrt P. Vector index is counterIJ
2176 
2177  //! Intermediate variable called Phi in Pitzer's paper. Vector index is
2178  //! counterIJ
2180 
2181  //! Derivative of m_Phi_IJ wrt T. Vector index is counterIJ
2183 
2184  //! Derivative of m_Phi_IJ wrt TT. Vector index is counterIJ
2186 
2187  //! Derivative of m_Phi_IJ wrt P. Vector index is counterIJ
2189 
2190  //! Intermediate variable called Phiprime in Pitzer's paper. Vector index is
2191  //! counterIJ
2193 
2194  //! Intermediate variable called PhiPhi in Pitzer's paper. Vector index is
2195  //! counterIJ
2197 
2198  //! Derivative of m_PhiPhi_IJ wrt T. Vector index is counterIJ
2200 
2201  //! Derivative of m_PhiPhi_IJ wrt TT. Vector index is counterIJ
2203 
2204  //! Derivative of m_PhiPhi_IJ wrt P. Vector index is counterIJ
2206 
2207  //! Intermediate variable called CMX in Pitzer's paper. Vector index is
2208  //! counterIJ
2210 
2211  //! Derivative of m_CMX_IJ wrt T. Vector index is counterIJ
2213 
2214  //! Derivative of m_CMX_IJ wrt TT. Vector index is counterIJ
2216 
2217  //! Derivative of m_CMX_IJ wrt P. Vector index is counterIJ
2219 
2220  //! Intermediate storage of the activity coefficient itself. Vector index is
2221  //! the species index
2223 
2224  //! Logarithm of the molal activity coefficients. Normally these are all
2225  //! one. However, stability schemes will change that
2227 
2228  //! value of the solute mole fraction that centers the cutoff polynomials
2229  //! for the cutoff =1 process;
2230  doublereal IMS_X_o_cutoff_;
2231 
2232  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
2233  doublereal IMS_cCut_;
2234 
2235  //! Parameter in the polyExp cutoff treatment
2236  /*!
2237  * This is the slope of the g function at the zero solvent point
2238  * Default value is 0.0
2239  */
2240  doublereal IMS_slopegCut_;
2241 
2242  //! @name Parameters in the polyExp cutoff treatment having to do with rate of exp decay
2243  //! @{
2244  doublereal IMS_dfCut_;
2245  doublereal IMS_efCut_;
2246  doublereal IMS_afCut_;
2247  doublereal IMS_bfCut_;
2248  doublereal IMS_dgCut_;
2249  doublereal IMS_egCut_;
2250  doublereal IMS_agCut_;
2251  doublereal IMS_bgCut_;
2252  //! @}
2253 
2254  //! value of the solvent mole fraction that centers the cutoff polynomials
2255  //! for the cutoff =1 process;
2256  doublereal MC_X_o_cutoff_;
2257 
2258  //! @name Parameters in the Molality Exp cutoff treatment
2259  //! @{
2260  doublereal MC_dpCut_;
2261  doublereal MC_epCut_;
2262  doublereal MC_apCut_;
2263  doublereal MC_bpCut_;
2264  doublereal MC_cpCut_;
2265  doublereal CROP_ln_gamma_o_min;
2266  doublereal CROP_ln_gamma_o_max;
2267  doublereal CROP_ln_gamma_k_min;
2268  doublereal CROP_ln_gamma_k_max;
2269 
2270  //! This is a boolean-type vector indicating whether
2271  //! a species's activity coefficient is in the cropped regime
2272  /*!
2273  * * 0 = Not in cropped regime
2274  * * 1 = In a transition regime where it is altered but there
2275  * still may be a temperature or pressure dependence
2276  * * 2 = In a cropped regime where there is no temperature
2277  * or pressure dependence
2278  */
2280  //! @}
2281 
2282  //! Initialize all of the species-dependent lengths in the object
2283  void initLengths();
2284 
2285  //! Apply the current phScale to a set of activity Coefficients or
2286  //! activities
2287  /*!
2288  * See the Eq3/6 Manual for a thorough discussion.
2289  *
2290  * @param acMolality input/Output vector containing the molality based
2291  * activity coefficients. length: m_kk.
2292  */
2293  virtual void applyphScale(doublereal* acMolality) const;
2294 
2295 private:
2296  /*
2297  * This function will be called to update the internally stored
2298  * natural logarithm of the molality activity coefficients
2299  */
2300  void s_update_lnMolalityActCoeff() const;
2301 
2302  //! This function calculates the temperature derivative of the
2303  //! natural logarithm of the molality activity coefficients.
2304  /*!
2305  * This function does all of the direct work. The solvent activity
2306  * coefficient is on the molality scale. It's derivative is too.
2307  */
2308  void s_update_dlnMolalityActCoeff_dT() const;
2309 
2310  /**
2311  * This function calculates the temperature second derivative of the natural
2312  * logarithm of the molality activity coefficients.
2313  */
2314  void s_update_d2lnMolalityActCoeff_dT2() const;
2315 
2316  /**
2317  * This function calculates the pressure derivative of the
2318  * natural logarithm of the molality activity coefficients.
2319  *
2320  * Assumes that the activity coefficients are current.
2321  */
2322  void s_update_dlnMolalityActCoeff_dP() const;
2323 
2324  //! This function will be called to update the internally stored
2325  //! natural logarithm of the molality activity coefficients
2326  /*
2327  * Normally they are all one. However, sometimes they are not,
2328  * due to stability schemes
2329  *
2330  * gamma_k_molar = gamma_k_molal / Xmol_solvent
2331  *
2332  * gamma_o_molar = gamma_o_molal
2333  */
2334  void s_updateIMS_lnMolalityActCoeff() const;
2335 
2336 private:
2337  //! Calculate the Pitzer portion of the activity coefficients.
2338  /**
2339  * This is the main routine in the whole module. It calculates the molality
2340  * based activity coefficients for the solutes, and the activity of water.
2341  */
2342  void s_updatePitzer_lnMolalityActCoeff() const;
2343 
2344  //! Calculates the temperature derivative of the natural logarithm of the
2345  //! molality activity coefficients.
2346  /*!
2347  * Public function makes sure that all dependent data is
2348  * up to date, before calling a private function
2349  */
2351 
2352  /**
2353  * This function calculates the temperature second derivative of the
2354  * natural logarithm of the molality activity coefficients.
2355  *
2356  * It is assumed that the Pitzer activity coefficient and first derivative
2357  * routine are called immediately preceding the call to this routine.
2358  */
2360 
2361  //! Calculates the Pressure derivative of the natural logarithm of the
2362  //! molality activity coefficients.
2363  /*!
2364  * It is assumed that the Pitzer activity coefficient and first derivative
2365  * routine are called immediately preceding the calling of this routine.
2366  */
2368 
2369  //! Calculates the Pitzer coefficients' dependence on the temperature.
2370  /*!
2371  * It will also calculate the temperature derivatives of the coefficients,
2372  * as they are important in the calculation of the latent heats and the heat
2373  * capacities of the mixtures.
2374  *
2375  * @param doDerivs If >= 1, then the routine will calculate the first
2376  * derivative. If >= 2, the routine will calculate the first
2377  * and second temperature derivative. default = 2
2378  */
2379  void s_updatePitzer_CoeffWRTemp(int doDerivs = 2) const;
2380 
2381  //! Calculate the lambda interactions.
2382  /*!
2383  * Calculate E-lambda terms for charge combinations of like sign, using
2384  * method of Pitzer (1975). This implementation is based on Bethke,
2385  * Appendix 2.
2386  *
2387  * @param is Ionic strength
2388  */
2389  void calc_lambdas(double is) const;
2390  mutable doublereal m_last_is;
2391 
2392  /**
2393  * Calculate etheta and etheta_prime
2394  *
2395  * This interaction accounts for the mixing effects of like-signed ions with
2396  * different charges. This interaction will be nonzero for species with the
2397  * same charge. this routine is not to be called for neutral species; it
2398  * core dumps or error exits.
2399  *
2400  * MEC implementation routine.
2401  *
2402  * @param z1 charge of the first molecule
2403  * @param z2 charge of the second molecule
2404  * @param etheta return pointer containing etheta
2405  * @param etheta_prime Return pointer containing etheta_prime.
2406  *
2407  * This routine uses the internal variables, elambda[] and elambda1[].
2408  */
2409  void calc_thetas(int z1, int z2,
2410  double* etheta, double* etheta_prime) const;
2411 
2412  //! Set up a counter variable for keeping track of symmetric binary
2413  //! interactions amongst the solute species.
2414  /*!
2415  * The purpose of this is to squeeze the ij parameters into a
2416  * compressed single counter.
2417  *
2418  * n = m_kk*i + j
2419  * m_Counter[n] = counter
2420  */
2421  void counterIJ_setup() const;
2422 
2423  //! Calculate the cropped molalities
2424  /*!
2425  * This is an internal routine that calculates values of m_molalitiesCropped
2426  * from m_molalities
2427  */
2428  void calcMolalitiesCropped() const;
2429 
2430  //! Process an XML node called "binarySaltParameters"
2431  /*!
2432  * This node contains all of the parameters necessary to describe the Pitzer
2433  * model for that particular binary salt. This function reads the XML file
2434  * and writes the coefficients it finds to an internal data structures.
2435  *
2436  * @param BinSalt reference to the XML_Node named binarySaltParameters
2437  * containing the anion - cation interaction
2438  */
2439  void readXMLBinarySalt(XML_Node& BinSalt);
2440 
2441  //! Process an XML node called "thetaAnion" or "thetaCation"
2442  /*!
2443  * This node contains all of the parameters necessary to describe the binary
2444  * interactions between two anions or two cations.
2445  *
2446  * @param BinSalt reference to the XML_Node named thetaAnion containing the
2447  * anion - anion interaction
2448  */
2449  void readXMLTheta(XML_Node& BinSalt);
2450 
2451  //! Process an XML node called "psiCommonAnion" or "psiCommonCation"
2452  /*!
2453  * This node contains all of the parameters necessary to describe
2454  * the ternary interactions between one anion and two cations or two anions
2455  * and one cation.
2456  *
2457  * @param BinSalt reference to the XML_Node named psiCommonAnion containing
2458  * the anion - cation1 - cation2 interaction
2459  */
2460  void readXMLPsi(XML_Node& BinSalt);
2461 
2462  //! Process an XML node called "lambdaNeutral"
2463  /*!
2464  * This node contains all of the parameters necessary to describe the binary
2465  * interactions between one neutral species and any other species (neutral
2466  * or otherwise) in the mechanism.
2467  *
2468  * @param BinSalt reference to the XML_Node named lambdaNeutral containing
2469  * multiple Neutral - species interactions
2470  */
2471  void readXMLLambdaNeutral(XML_Node& BinSalt);
2472 
2473  //! Process an XML node called "MunnnNeutral"
2474  /*!
2475  * This node contains all of the parameters necessary to describe
2476  * the self-ternary interactions for one neutral species.
2477  *
2478  * @param BinSalt reference to the XML_Node named Munnn containing the
2479  * self-ternary interaction
2480  */
2481  void readXMLMunnnNeutral(XML_Node& BinSalt);
2482 
2483  //! Process an XML node called "zetaCation"
2484  /*!
2485  * This node contains all of the parameters necessary to describe
2486  * the ternary interactions between one neutral, one cation, and one anion.
2487  *
2488  * @param BinSalt reference to the XML_Node named psiCommonCation
2489  * containing the neutral - cation - anion interaction
2490  */
2491  void readXMLZetaCation(const XML_Node& BinSalt);
2492 
2493  //! Precalculate the IMS Cutoff parameters for typeCutoff = 2
2494  void calcIMSCutoffParams_();
2495 
2496  //! Calculate molality cut-off parameters
2497  void calcMCCutoffParams_();
2498 };
2499 
2500 }
2501 
2502 #endif
vector_fp m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1927
vector_fp m_BMX_IJ_L
Derivative of BMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2143
vector_fp m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2205
void calcMolalitiesCropped() const
Calculate the cropped molalities.
Definition: HMWSoln.cpp:1172
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: HMWSoln.cpp:312
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
Definition: HMWSoln.cpp:4109
vector_fp m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1894
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
Definition: HMWSoln.cpp:2399
vector_fp m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1921
void readXMLBinarySalt(XML_Node &BinSalt)
Process an XML node called "binarySaltParameters".
Definition: HMWSoln.cpp:1336
vector_fp m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2218
vector_fp m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
Definition: HMWSoln.h:2125
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
Definition: HMWSoln.cpp:210
vector_fp m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1850
vector_fp m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
Definition: HMWSoln.h:2085
doublereal MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
Definition: HMWSoln.h:2256
virtual void applyphScale(doublereal *acMolality) const
Apply the current phScale to a set of activity Coefficients or activities.
Definition: HMWSoln.cpp:4095
vector_fp IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
Definition: HMWSoln.h:2226
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2014
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition: HMWSoln.h:1726
void readXMLLambdaNeutral(XML_Node &BinSalt)
Process an XML node called "lambdaNeutral".
Definition: HMWSoln.cpp:1471
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
Definition: HMWSoln.cpp:198
vector_fp m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1829
doublereal s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
Definition: HMWSoln.cpp:4191
vector_fp m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2032
vector_fp m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2093
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer&#39;s papers.
Definition: HMWSoln.h:1816
vector_fp m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1856
vector_fp m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1847
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
Definition: HMWSoln.h:1999
void readXMLMunnnNeutral(XML_Node &BinSalt)
Process an XML node called "MunnnNeutral".
Definition: HMWSoln.cpp:1501
vector_fp m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer&#39;s paper.
Definition: HMWSoln.h:2196
vector_fp m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2182
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
Definition: HMWSoln.cpp:351
vector_fp m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2073
virtual void setMolarDensity(const doublereal conc)
Set the internally stored molar density (kmol/m^3) for the phase.
Definition: HMWSoln.cpp:252
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
Definition: HMWSoln.cpp:1313
vector_fp m_BMX_IJ_P
Derivative of BMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2149
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition: HMWSoln.cpp:4008
vector_fp m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1802
double elambda[17]
This is elambda, MEC.
Definition: HMWSoln.h:2113
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
vector_fp m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2089
vector_fp m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2052
doublereal s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale...
Definition: HMWSoln.cpp:4177
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
Definition: HMWSoln.h:1993
virtual void setDensity(const doublereal rho)
Set the internally stored density (kg/m^3) of the phase.
Definition: HMWSoln.cpp:243
vector_fp m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2162
void calcMCCutoffParams_()
Calculate molality cut-off parameters.
Definition: HMWSoln.cpp:1614
vector_fp m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1897
void calcIMSCutoffParams_()
Precalculate the IMS Cutoff parameters for typeCutoff = 2.
Definition: HMWSoln.cpp:1561
vector_fp m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
Definition: HMWSoln.h:1956
void readXMLZetaCation(const XML_Node &BinSalt)
Process an XML node called "zetaCation".
Definition: HMWSoln.cpp:1526
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Initialize the phase parameters from an XML file.
Definition: HMWSoln.cpp:746
vector_fp m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2069
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:2946
vector_fp m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1918
vector_fp m_hfunc_IJ
hfunc, was called gprime in Pitzer&#39;s paper.
Definition: HMWSoln.h:2132
double elambda1[17]
This is elambda1, MEC.
Definition: HMWSoln.h:2116
vector_fp m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1832
vector_fp m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2159
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer&#39;s definition of A_V.
Definition: HMWSoln.cpp:944
vector_fp m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2202
doublereal IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
Definition: HMWSoln.h:2233
vector_fp m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
Definition: HMWSoln.h:1960
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
Definition: HMWSoln.h:1722
vector_fp m_Phi_IJ
Intermediate variable called Phi in Pitzer&#39;s paper.
Definition: HMWSoln.h:2179
vector_fp m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2166
void setA_Debye(double A)
Set the A_Debye parameter.
Definition: HMWSoln.cpp:641
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
Definition: HMWSoln.h:1906
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
Definition: HMWSoln.cpp:386
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure...
Definition: HMWSoln.cpp:4154
vector_fp m_Beta1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1823
vector_fp m_tmpV
vector of size m_kk, used as a temporary holding area.
Definition: HMWSoln.h:1790
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:31
Header for intermediate ThermoPhase object for phases which employ molality based activity coefficien...
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
Definition: HMWSoln.cpp:399
HMWSoln()
Default Constructor.
Definition: HMWSoln.cpp:30
vector_fp m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2062
Header file for class Cantera::Array2D.
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
vector_fp m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2188
PDSS * m_waterSS
Water standard state calculator.
Definition: HMWSoln.h:1784
void readXMLTheta(XML_Node &BinSalt)
Process an XML node called "thetaAnion" or "thetaCation".
Definition: HMWSoln.cpp:1382
vector_fp m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer&#39;s paper.
Definition: HMWSoln.h:2192
vector_fp m_h2func_IJ
hfunc2, was called gprime in Pitzer&#39;s paper.
Definition: HMWSoln.h:2136
void printCoeffs() const
Print out all of the input Pitzer coefficients.
Definition: HMWSoln.cpp:4051
vector_fp m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
Definition: HMWSoln.h:1964
vector_fp m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2212
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
Definition: HMWSoln.h:2102
vector_fp m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2077
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1980
doublereal s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
Definition: HMWSoln.cpp:4169
std::unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition: HMWSoln.h:1787
virtual doublereal relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
Definition: HMWSoln.cpp:146
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
Definition: HMWSoln.cpp:3982
virtual void initThermo()
Definition: HMWSoln.cpp:660
vector_fp m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2022
vector_fp m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
Definition: HMWSoln.h:2099
vector_fp m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1853
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients...
Definition: HMWSoln.cpp:2427
vector_int CROP_speciesCropped_
This is a boolean-type vector indicating whether a species&#39;s activity coefficient is in the cropped r...
Definition: HMWSoln.h:2279
vector_fp m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2215
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer&#39;s papers.
Definition: HMWSoln.h:1866
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
Definition: HMWSoln.h:1996
void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
Definition: HMWSoln.cpp:299
doublereal IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
Definition: HMWSoln.h:2240
vector_fp m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1805
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients&#39; dependence on the temperature.
Definition: HMWSoln.cpp:1650
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
Definition: HMWSoln.h:1746
doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale...
Definition: HMWSoln.cpp:4184
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
Definition: HMWSoln.cpp:140
vector_int m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species...
Definition: HMWSoln.h:2110
vector_fp m_CMX_IJ
Intermediate variable called CMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2209
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
Definition: HMWSoln.cpp:332
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
Definition: HMWSoln.h:1729
vector_fp m_CphiMX_ij_L
Derivative of Cphi_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1891
void readXMLPsi(XML_Node &BinSalt)
Process an XML node called "psiCommonAnion" or "psiCommonCation".
Definition: HMWSoln.cpp:1422
vector_fp m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2169
vector_fp m_g2func_IJ
This is the value of g2(x2) in Pitzer&#39;s papers. Vector index is counterIJ.
Definition: HMWSoln.h:2128
virtual double d2A_DebyedT2_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the 2nd derivative of the Debye Huckel constant with respect to temperature as a function of...
Definition: HMWSoln.cpp:967
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
Definition: HMWSoln.h:1990
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:2918
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
Definition: HMWSoln.cpp:216
vector_fp m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2185
virtual doublereal pressure() const
Returns the current pressure of the phase.
virtual std::string type() const
String indicating the thermodynamic model implemented.
Definition: HMWSoln.h:1070
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized activity concentrations.
Definition: HMWSoln.cpp:260
vector_fp m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1881
vector_fp m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1799
vector_fp m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2156
vector_fp m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2146
double m_A_Debye
A_Debye: this expression appears on the top of the ln actCoeff term in the general Debye-Huckel expre...
Definition: HMWSoln.h:1778
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
Definition: HMWSoln.cpp:273
void initLengths()
Initialize all of the species-dependent lengths in the object.
Definition: HMWSoln.cpp:995
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
Virtual base class for a species with a pressure dependent standard state.
Definition: PDSS.h:165
virtual doublereal relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
Definition: HMWSoln.cpp:158
void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
Definition: HMWSoln.cpp:228
vector_fp m_gamma_tmp
Intermediate storage of the activity coefficient itself.
Definition: HMWSoln.h:2222
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
Definition: HMWSoln.cpp:3429
virtual void getActivities(doublereal *ac) const
Get the array of non-dimensional activities at the current solution temperature, pressure, and solution concentration.
Definition: HMWSoln.cpp:283
doublereal IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
Definition: HMWSoln.h:2230
vector_fp m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1952
Class HMWSoln represents a dilute or concentrated liquid electrolyte phase which obeys the Pitzer for...
Definition: HMWSoln.h:1038
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients...
Definition: HMWSoln.cpp:3453
vector_fp m_BMX_IJ
Intermediate variable called BMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2140
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature...
Definition: HMWSoln.cpp:4124
virtual double dA_DebyedP_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to pressure, as a function of tempe...
Definition: HMWSoln.cpp:901
virtual double dA_DebyedT_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to temperature as a function of tem...
Definition: HMWSoln.cpp:877
void s_updateScaling_pHScaling_dT2() const
Apply the current phScale to a set of 2nd derivatives of the activity Coefficients wrt temperature...
Definition: HMWSoln.cpp:4139
vector_fp m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2153
vector_fp m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
Definition: HMWSoln.h:2081
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
vector_fp m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1888
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model...
Definition: HMWSoln.h:1717
void calc_lambdas(double is) const
Calculate the lambda interactions.
Definition: HMWSoln.cpp:3938
vector_fp m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1826
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer&#39;s definition of A_J.
Definition: HMWSoln.cpp:955
vector_fp m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1924
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1939
virtual doublereal satPressure(doublereal T)
Get the saturation pressure for a given temperature.
Definition: HMWSoln.cpp:419
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
Definition: HMWSoln.h:2055
vector_fp m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2042
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
Definition: HMWSoln.cpp:1896
vector_fp m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2175
virtual double A_Debye_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the Debye Huckel constant as a function of temperature and pressure.
Definition: HMWSoln.cpp:845
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer&#39;s papers.
Definition: HMWSoln.h:1840
vector_fp m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2199
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer&#39;s definition of A_L.
Definition: HMWSoln.cpp:933
vector_fp m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1808
vector_fp m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2172
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Definition: HMWSoln.cpp:204