Cantera  2.0
DebyeHuckel.h
Go to the documentation of this file.
1 /**
2  * @file DebyeHuckel.h
3  * Headers for the %DebyeHuckel ThermoPhase object, which models dilute
4  * electrolyte solutions
5  * (see \ref thermoprops and \link Cantera::DebyeHuckel DebyeHuckel \endlink) .
6  *
7  * Class %DebyeHuckel represents a dilute liquid electrolyte phase which
8  * obeys the Debye Huckel formulation for nonideality.
9  */
10 
11 /*
12  * Copyright (2006) Sandia Corporation. Under the terms of
13  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
14  * U.S. Government retains certain rights in this software.
15  */
16 #ifndef CT_DEBYEHUCKEL_H
17 #define CT_DEBYEHUCKEL_H
18 
19 #include "MolalityVPSSTP.h"
20 #include "electrolytes.h"
21 #include "cantera/base/Array.h"
22 
23 namespace Cantera
24 {
25 
26 /*!
27  * @name Formats for the Activity Coefficients
28  *
29  * These are possible formats for the molality-based activity coefficients.
30  */
31 //@{
32 #define DHFORM_DILUTE_LIMIT 0
33 #define DHFORM_BDOT_AK 1
34 #define DHFORM_BDOT_ACOMMON 2
35 #define DHFORM_BETAIJ 3
36 #define DHFORM_PITZER_BETAIJ 4
37 //@}
38 /*
39  * @name Acceptable ways to calculate the value of A_Debye
40  */
41 //@{
42 #define A_DEBYE_CONST 0
43 #define A_DEBYE_WATER 1
44 //@}
45 
46 class WaterProps;
47 class PDSS_Water;
48 
49 /**
50  * @ingroup thermoprops
51  *
52  * Class %DebyeHuckel represents a dilute liquid electrolyte phase which
53  * obeys the Debye Huckel formulation for nonideality.
54  *
55  * The concentrations of the ionic species are assumed to obey the electroneutrality
56  * condition.
57  *
58  * <HR>
59  * <H2> Specification of Species Standard %State Properties </H2>
60  * <HR>
61  *
62  * The standard states are on the unit molality basis. Therefore, in the
63  * documentation below, the normal \f$ o \f$ superscript is replaced with
64  * the \f$ \triangle \f$ symbol. The reference state symbol is now
65  * \f$ \triangle, ref \f$.
66  *
67  *
68  * It is assumed that the reference state thermodynamics may be
69  * obtained by a pointer to a populated species thermodynamic property
70  * manager class (see ThermoPhase::m_spthermo). How to relate pressure
71  * changes to the reference state thermodynamics is resolved at this level.
72  *
73  * For an incompressible,
74  * stoichiometric substance, the molar internal energy is
75  * independent of pressure. Since the thermodynamic properties
76  * are specified by giving the standard-state enthalpy, the
77  * term \f$ P_0 \hat v\f$ is subtracted from the specified molar
78  * enthalpy to compute the molar internal energy. The entropy is
79  * assumed to be independent of the pressure.
80  *
81  * The enthalpy function is given by the following relation.
82  *
83  * \f[
84  * \raggedright h^\triangle_k(T,P) = h^{\triangle,ref}_k(T)
85  * + \tilde v \left( P - P_{ref} \right)
86  * \f]
87  *
88  * For an incompressible,
89  * stoichiometric substance, the molar internal energy is
90  * independent of pressure. Since the thermodynamic properties
91  * are specified by giving the standard-state enthalpy, the
92  * term \f$ P_{ref} \tilde v\f$ is subtracted from the specified reference molar
93  * enthalpy to compute the molar internal energy.
94  *
95  * \f[
96  * u^\triangle_k(T,P) = h^{\triangle,ref}_k(T) - P_{ref} \tilde v
97  * \f]
98  *
99  * The standard state heat capacity and entropy are independent
100  * of pressure. The standard state gibbs free energy is obtained
101  * from the enthalpy and entropy functions.
102  *
103  * The vector Phase::m_speciesSize[] is used to hold the
104  * base values of species sizes. These are defined as the
105  * molar volumes of species at infinite dilution at 300 K and 1 atm
106  * of water. m_speciesSize are calculated during the initialization of the
107  * %DebyeHuckel object and are then not touched.
108  *
109  * The current model assumes that an incompressible molar volume for
110  * all solutes. The molar volume for the water solvent, however,
111  * is obtained from a pure water equation of state, waterSS.
112  * Therefore, the water standard state varies with both T and P.
113  * It is an error to request standard state water properties at a T and P
114  * where the water phase is not a stable phase, i.e., beyond its
115  * spinodal curve.
116  *
117  * <HR>
118  * <H2> Specification of Solution Thermodynamic Properties </H2>
119  * <HR>
120  *
121  * Chemical potentials
122  * of the solutes, \f$ \mu_k \f$, and the solvent, \f$ \mu_o \f$, which are based
123  * on the molality form, have the following general format:
124  *
125  * \f[
126  * \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} \frac{m_k}{m^\triangle})
127  * \f]
128  * \f[
129  * \mu_o = \mu^o_o(T,P) + RT ln(a_o)
130  * \f]
131  *
132  * where \f$ \gamma_k^{\triangle} \f$ is the molality based activity coefficient for species
133  * \f$k\f$.
134  *
135  * Individual activity coefficients of ions can not be independently measured. Instead,
136  * only binary pairs forming electroneutral solutions can be measured.
137 
138  *
139  * <H3> Ionic Strength </H3>
140  *
141  * Most of the parameterizations within the model use the ionic strength
142  * as a key variable. The ionic strength, \f$ I\f$ is defined as follows
143  *
144  * \f[
145  * I = \frac{1}{2} \sum_k{m_k z_k^2}
146  * \f]
147  *
148  * \f$ m_k \f$ is the molality of the kth species. \f$ z_k \f$ is the charge
149  * of the kth species. Note, the ionic strength is a defined units quantity.
150  * The molality has defined units of gmol kg-1, and therefore the ionic
151  * strength has units of sqrt( gmol kg-1).
152  *
153  * In some instances, from some authors, a different
154  * formulation is used for the ionic strength in the equations below. The different
155  * formulation is due to the possibility of the existence of weak acids and how
156  * association wrt to the weak acid equilibrium relation affects the calculation
157  * of the activity coefficients via the assumed value of the ionic strength.
158  *
159  * If we are to assume that the association reaction doesn't have an effect
160  * on the ionic strength, then we will want to consider the associated weak
161  * acid as in effect being fully dissociated, when we calculate an effective
162  * value for the ionic strength. We will call this calculated value, the
163  * stoichiometric ionic strength, \f$ I_s \f$, putting a subscript s to denote
164  * it from the more straightforward calculation of \f$ I \f$.
165  *
166  * \f[
167  * I_s = \frac{1}{2} \sum_k{m_k^s z_k^2}
168  * \f]
169  *
170  * Here, \f$ m_k^s \f$ is the value of the molalities calculated assuming that
171  * all weak acid-base pairs are in their fully dissociated states. This calculation may
172  * be simplified by considering that the weakly associated acid may be made up of two
173  * charged species, k1 and k2, each with their own charges, obeying the following relationship:
174  *
175  * \f[
176  * z_k = z_{k1} + z_{k2}
177  * \f]
178  * Then, we may only need to specify one charge value, say, \f$ z_{k1}\f$,
179  * the cation charge number,
180  * in order to get both numbers, since we have already specified \f$ z_k \f$
181  * in the definition of original species.
182  * Then, the stoichiometric ionic strength may be calculated via the following formula.
183  *
184  * \f[
185  * I_s = \frac{1}{2} \left(\sum_{k,ions}{m_k z_k^2}+
186  * \sum_{k,weak_assoc}(m_k z_{k1}^2 + m_k z_{k2}^2) \right)
187  * \f]
188  *
189  * The specification of which species are weakly associated acids is made in the input
190  * file via the
191  * <TT> stoichIsMods </TT> XML block, where the charge for k1 is also specified.
192  * An example is given below:
193  *
194  * @code
195  * <stoichIsMods>
196  * NaCl(aq):-1.0
197  * </stoichIsMods>
198  * @endcode
199  *
200  * Because we need the concept of a weakly associated acid in order to calculated
201  * \f$ I_s \f$ we need to
202  * catalog all species in the phase. This is done using the following categories:
203  *
204  * - <B>cEST_solvent</B> : Solvent species (neutral)
205  * - <B>cEST_chargedSpecies</B> Charged species (charged)
206  * - <B>cEST_weakAcidAssociated</B> Species which can break apart into charged species.
207  * It may or may not be charged. These may or
208  * may not be be included in the
209  * species solution vector.
210  * - <B>cEST_strongAcidAssociated</B> Species which always breaks apart into charged species.
211  * It may or may not be charged. Normally, these aren't included
212  * in the speciation vector.
213  * - <B>cEST_polarNeutral </B> Polar neutral species
214  * - <B>cEST_nonpolarNeutral</B> Non polar neutral species
215  *
216  * Polar and non-polar neutral species are differentiated, because some additions
217  * to the activity
218  * coefficient expressions distinguish between these two types of solutes. This is the so-called
219  * salt-out effect.
220  *
221  * The type of species is specified in the <TT>electrolyteSpeciesType</TT> XML block.
222  * Note, this is not
223  * considered a part of the specification of the standard state for the species,
224  * at this time. Therefore,
225  * this information is put under the <TT>activityCoefficient</TT> XML block. An example
226  * is given below
227  *
228  * @code
229  * <electrolyteSpeciesType>
230  * H2L(L):solvent
231  * H+:chargedSpecies
232  * NaOH(aq):weakAcidAssociated
233  * NaCl(aq):strongAcidAssociated
234  * NH3(aq):polarNeutral
235  * O2(aq):nonpolarNeutral
236  * </electrolyteSpeciesType>
237  * @endcode
238  *
239  * Much of the species electrolyte type information is inferred from other information in the
240  * input file. For example, as species which is charged is given the "chargedSpecies" default
241  * category. A neutral solute species is put into the "nonpolarNeutral" category by default.
242  *
243  * The specification of solute activity coefficients depends on the model
244  * assumed for the Debye-Huckel term. The model is set by the
245  * internal parameter #m_formDH. We will now describe each category in its own section.
246  *
247  *
248  * <H3> Debye-Huckel Dilute Limit </H3>
249  *
250  * DHFORM_DILUTE_LIMIT = 0
251  *
252  * This form assumes a dilute limit to DH, and is mainly
253  * for informational purposes:
254  * \f[
255  * \ln(\gamma_k^\triangle) = - z_k^2 A_{Debye} \sqrt{I}
256  * \f]
257  * where \f$ I\f$ is the ionic strength
258  * \f[
259  * I = \frac{1}{2} \sum_k{m_k z_k^2}
260  * \f]
261  *
262  * The activity for the solvent water,\f$ a_o \f$, is not independent and must be
263  * determined from the Gibbs-Duhem relation.
264  *
265  * \f[
266  * \ln(a_o) = \frac{X_o - 1.0}{X_o} + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{3/2}
267  * \f]
268  *
269  *
270  * <H3> Bdot Formulation </H3>
271  *
272  * DHFORM_BDOT_AK = 1
273  *
274  * This form assumes Bethke's format for the Debye Huckel activity coefficient:
275  *
276  * \f[
277  * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye} \sqrt{I}}{ 1 + B_{Debye} a_k \sqrt{I}}
278  * + \log(10) B^{dot}_k I
279  * \f]
280  *
281  * Note, this particular form where \f$ a_k \f$ can differ in
282  * multielectrolyte
283  * solutions has problems with respect to a Gibbs-Duhem analysis. However,
284  * we include it here because there is a lot of data fit to it.
285  *
286  * The activity for the solvent water,\f$ a_o \f$, is not independent and must be
287  * determined from the Gibbs-Duhem relation. Here, we use:
288  *
289  * \f[
290  * \ln(a_o) = \frac{X_o - 1.0}{X_o}
291  * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{1/2}
292  * \left[ \sum_k{\frac{1}{2} m_k z_k^2 \sigma( B_{Debye} a_k \sqrt{I} ) } \right]
293  * - \frac{\log(10)}{2} \tilde{M}_o I \sum_k{ B^{dot}_k m_k}
294  * \f]
295  * where
296  * \f[
297  * \sigma (y) = \frac{3}{y^3} \left[ (1+y) - 2 \ln(1 + y) - \frac{1}{1+y} \right]
298  * \f]
299  *
300  * Additionally, Helgeson's formulation for the water activity is offered as an
301  * alternative.
302  *
303  *
304  * <H3> Bdot Formulation with Uniform Size Parameter in the Denominator </H3>
305  *
306  * DHFORM_BDOT_AUNIFORM = 2
307  *
308  * This form assumes Bethke's format for the Debye-Huckel activity coefficient
309  *
310  * \f[
311  * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye} \sqrt{I}}{ 1 + B_{Debye} a \sqrt{I}}
312  * + \log(10) B^{dot}_k I
313  * \f]
314  *
315  * The value of a is determined at the beginning of the
316  * calculation, and not changed.
317  *
318  * \f[
319  * \ln(a_o) = \frac{X_o - 1.0}{X_o}
320  * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{3/2} \sigma( B_{Debye} a \sqrt{I} )
321  * - \frac{\log(10)}{2} \tilde{M}_o I \sum_k{ B^{dot}_k m_k}
322  * \f]
323  *
324  *
325  * <H3> Beta_IJ formulation </H3>
326  *
327  * DHFORM_BETAIJ = 3
328  *
329  * This form assumes a linear expansion in a virial coefficient form
330  * It is used extensively in the book by Newmann, "Electrochemistry Systems",
331  * and is the beginning of
332  * more complex treatments for stronger electrolytes, fom Pitzer
333  * and from Harvey, Moller, and Weire.
334  *
335  * \f[
336  * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye} \sqrt{I}}{ 1 + B_{Debye} a \sqrt{I}}
337  * + 2 \sum_j \beta_{j,k} m_j
338  * \f]
339  *
340  * In the current treatment the binary interaction coefficients, \f$ \beta_{j,k}\f$, are
341  * independent of temperature and pressure.
342  *
343  * \f[
344  * \ln(a_o) = \frac{X_o - 1.0}{X_o}
345  * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{3/2} \sigma( B_{Debye} a \sqrt{I} )
346  * - \tilde{M}_o \sum_j \sum_k \beta_{j,k} m_j m_k
347  * \f]
348  *
349  * In this formulation the ionic radius, \f$ a \f$, is a constant. This must be supplied to the
350  * model, in an <DFN> ionicRadius </DFN> XML block.
351  *
352  * The \f$ \beta_{j,k} \f$ parameters are binary interaction parameters. They are supplied to
353  * the object in an <TT> DHBetaMatrix </TT> XML block. There are in principle \f$ N (N-1) /2 \f$
354  * different, symmetric interaction parameters, where \f$ N \f$ are the number of solute species in the
355  * mechanism.
356  * An example is given below.
357  *
358  * An example <TT> activityCoefficients </TT> XML block for this formulation is supplied below
359  *
360  * @code
361  * <activityCoefficients model="Beta_ij">
362  * <!-- A_Debye units = sqrt(kg/gmol) -->
363  * <A_Debye> 1.172576 </A_Debye>
364  * <!-- B_Debye units = sqrt(kg/gmol)/m -->
365  * <B_Debye> 3.28640E9 </B_Debye>
366  * <ionicRadius default="3.042843" units="Angstroms">
367  * </ionicRadius>
368  * <DHBetaMatrix>
369  * H+:Cl-:0.27
370  * Na+:Cl-:0.15
371  * Na+:OH-:0.06
372  * </DHBetaMatrix>
373  * <stoichIsMods>
374  * NaCl(aq):-1.0
375  * </stoichIsMods>
376  * <electrolyteSpeciesType>
377  * H+:chargedSpecies
378  * NaCl(aq):weakAcidAssociated
379  * </electrolyteSpeciesType>
380  * </activityCoefficients>
381  * @endcode
382  *
383  * <H3> Pitzer Beta_IJ formulation </H3>
384  *
385  * DHFORM_PITZER_BETAIJ = 4
386  *
387  * This form assumes an activity coefficient formulation consistent
388  * with a truncated form of Pitzer's formulation. Pitzer's formulation is equivalent
389  * to the formulations above in the dilute limit, where rigorous theory may be applied.
390  *
391  * \f[
392  * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye}}{3} \frac{\sqrt{I}}{ 1 + B_{Debye} a \sqrt{I}}
393  * -2 z_k^2 \frac{A_{Debye}}{3} \frac{\ln(1 + B_{Debye} a \sqrt{I})}{ B_{Debye} a}
394  * + 2 \sum_j \beta_{j,k} m_j
395  * \f]
396  *
397  *
398  * \f[
399  * \ln(a_o) = \frac{X_o - 1.0}{X_o}
400  * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} \frac{(I)^{3/2} }{1 + B_{Debye} a \sqrt{I} }
401  * - \tilde{M}_o \sum_j \sum_k \beta_{j,k} m_j m_k
402  * \f]
403  *
404  * <H3> Specification of the Debye Huckel Constants </H3>
405  *
406  * In the equations above, the formulas for \f$ A_{Debye} \f$ and \f$ B_{Debye} \f$
407  * are needed. The %DebyeHuckel object uses two methods for specifying these quantities.
408  * The default method is to assume that \f$ A_{Debye} \f$ is a constant, given
409  * in the initialization process, and stored in the
410  * member double, m_A_Debye. Optionally, a full water treatment may be employed that makes
411  * \f$ A_{Debye} \f$ a full function of <I>T</I> and <I>P</I>.
412  *
413  * \f[
414  * A_{Debye} = \frac{F e B_{Debye}}{8 \pi \epsilon R T} {\left( C_o \tilde{M}_o \right)}^{1/2}
415  * \f]
416  * where
417  *
418  * \f[
419  * B_{Debye} = \frac{F} {{(\frac{\epsilon R T}{2})}^{1/2}}
420  * \f]
421  * Therefore:
422  * \f[
423  * A_{Debye} = \frac{1}{8 \pi}
424  * {\left(\frac{2 N_a \rho_o}{1000}\right)}^{1/2}
425  * {\left(\frac{N_a e^2}{\epsilon R T }\right)}^{3/2}
426  * \f]
427  *
428  * Units = sqrt(kg/gmol)
429  *
430  * where
431  * - \f$ N_a \f$ is Avogadro's number
432  * - \f$ \rho_w \f$ is the density of water
433  * - \f$ e \f$ is the electronic charge
434  * - \f$ \epsilon = K \epsilon_o \f$ is the permittivity of water
435  * where \f$ K \f$ is the dielectric constant of water,
436  * and \f$ \epsilon_o \f$ is the permittivity of free space.
437  * - \f$ \rho_o \f$ is the density of the solvent in its standard state.
438  *
439  * Nominal value at 298 K and 1 atm = 1.172576 (kg/gmol)<SUP>1/2</SUP>
440  * based on:
441  * - \f$ \epsilon / \epsilon_0 \f$ = 78.54
442  * (water at 25C)
443  * - T = 298.15 K
444  * - B_Debye = 3.28640E9 (kg/gmol)<SUP>1/2</SUP> m<SUP>-1</SUP>
445  *
446  * An example of a fixed value implementation is given below.
447  * @code
448  * <activityCoefficients model="Beta_ij">
449  * <!-- A_Debye units = sqrt(kg/gmol) -->
450  * <A_Debye> 1.172576 </A_Debye>
451  * <!-- B_Debye units = sqrt(kg/gmol)/m -->
452  * <B_Debye> 3.28640E9 </B_Debye>
453  * </activityCoefficients>
454  * @endcode
455  *
456  * An example of a variable value implementation is given below.
457  *
458  * @code
459  * <activityCoefficients model="Beta_ij">
460  * <A_Debye model="water" />
461  * <!-- B_Debye units = sqrt(kg/gmol)/m -->
462  * <B_Debye> 3.28640E9 </B_Debye>
463  * </activityCoefficients>
464  * @endcode
465  *
466  * Currently, \f$ B_{Debye} \f$ is a constant in the model, specified either by a default
467  * water value, or through the input file. This may have to be looked at, in the future.
468  *
469  * <HR>
470  * <H2> %Application within %Kinetics Managers </H2>
471  * <HR>
472  *
473  * For the time being, we have set the standard concentration for all species in
474  * this phase equal to the default concentration of the solvent at 298 K and 1 atm.
475  * This means that the
476  * kinetics operator essentially works on an activities basis, with units specified
477  * as if it were on a concentration basis.
478  *
479  * For example, a bulk-phase binary reaction between liquid species j and k, producing
480  * a new liquid species l would have the
481  * following equation for its rate of progress variable, \f$ R^1 \f$, which has
482  * units of kmol m-3 s-1.
483  *
484  * \f[
485  * R^1 = k^1 C_j^a C_k^a = k^1 (C_o a_j) (C_o a_k)
486  * \f]
487  * where
488  * \f[
489  * C_j^a = C_o a_j \quad and \quad C_k^a = C_o a_k
490  * \f]
491  *
492  * \f$ C_j^a \f$ is the activity concentration of species j, and
493  * \f$ C_k^a \f$ is the activity concentration of species k. \f$ C_o \f$
494  * is the concentration of water at 298 K and 1 atm. \f$ a_j \f$ is
495  * the activity of species j at the current temperature and pressure
496  * and concentration of the liquid phase. \f$k^1 \f$ has units of m3 kmol-1 s-1.
497  *
498  * The reverse rate constant can then be obtained from the law of microscopic reversibility
499  * and the equilibrium expression for the system.
500  *
501  * \f[
502  * \frac{a_j a_k}{ a_l} = K^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} )
503  * \f]
504  *
505  * \f$ K^{o,1} \f$ is the dimensionless form of the equilibrium constant.
506  *
507  * \f[
508  * R^{-1} = k^{-1} C_l^a = k^{-1} (C_o a_l)
509  * \f]
510  *
511  * where
512  *
513  * \f[
514  * k^{-1} = k^1 K^{o,1} C_o
515  * \f]
516  *
517  * \f$k^{-1} \f$ has units of s-1.
518  *
519  * Note, this treatment may be modified in the future, as events dictate.
520  *
521  * <HR>
522  * <H2> Instantiation of the Class </H2>
523  * <HR>
524  *
525  * The constructor for this phase is NOT located in the default ThermoFactory
526  * for %Cantera. However, a new %DebyeHuckel object may be created by
527  * the following code snippets:
528  *
529  * @code
530  * DebyeHuckel *DH = new DebyeHuckel("DH_NaCl.xml", "NaCl_electrolyte");
531  * @endcode
532  *
533  * or
534  *
535  * @code
536  * char iFile[80], file_ID[80];
537  * strcpy(iFile, "DH_NaCl.xml");
538  * sprintf(file_ID,"%s#NaCl_electrolyte", iFile);
539  * XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
540  * DebyeHuckel *dh = new DebyeHuckel(*xm);
541  * @endcode
542  *
543  * or by the following call to importPhase():
544  *
545  * @code
546  * char iFile[80], file_ID[80];
547  * strcpy(iFile, "DH_NaCl.xml");
548  * sprintf(file_ID,"%s#NaCl_electrolyte", iFile);
549  * XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
550  * DebyeHuckel dhphase;
551  * importPhase(*xm, &dhphase);
552  * @endcode
553  *
554  * <HR>
555  * <H2> XML Example </H2>
556  * <HR>
557  *
558  * The phase model name for this is called StoichSubstance. It must be supplied
559  * as the model attribute of the thermo XML element entry.
560  * Within the phase XML block,
561  * the density of the phase must be specified. An example of an XML file
562  * this phase is given below.
563  *
564  * @verbatim
565  <phase id="NaCl_electrolyte" dim="3">
566  <speciesArray datasrc="#species_waterSolution">
567  H2O(L) Na+ Cl- H+ OH- NaCl(aq) NaOH(aq)
568  </speciesArray>
569  <state>
570  <temperature units="K"> 300 </temperature>
571  <pressure units="Pa">101325.0</pressure>
572  <soluteMolalities>
573  Na+:3.0
574  Cl-:3.0
575  H+:1.0499E-8
576  OH-:1.3765E-6
577  NaCl(aq):0.98492
578  NaOH(aq):3.8836E-6
579  </soluteMolalities>
580  </state>
581  <!-- thermo model identifies the inherited class
582  from ThermoPhase that will handle the thermodynamics.
583  -->
584  <thermo model="DebyeHuckel">
585  <standardConc model="solvent_volume" />
586  <activityCoefficients model="Beta_ij">
587  <!-- A_Debye units = sqrt(kg/gmol) -->
588  <A_Debye> 1.172576 </A_Debye>
589  <!-- B_Debye units = sqrt(kg/gmol)/m -->
590  <B_Debye> 3.28640E9 </B_Debye>
591  <ionicRadius default="3.042843" units="Angstroms">
592  </ionicRadius>
593  <DHBetaMatrix>
594  H+:Cl-:0.27
595  Na+:Cl-:0.15
596  Na+:OH-:0.06
597  </DHBetaMatrix>
598  <stoichIsMods>
599  NaCl(aq):-1.0
600  </stoichIsMods>
601  <electrolyteSpeciesType>
602  H+:chargedSpecies
603  NaCl(aq):weakAcidAssociated
604  </electrolyteSpeciesType>
605  </activityCoefficients>
606  <solvent> H2O(L) </solvent>
607  </thermo>
608  <elementArray datasrc="elements.xml"> O H Na Cl </elementArray>
609 </phase>
610 @endverbatim
611  *
612  *
613  */
614 class DebyeHuckel : public MolalityVPSSTP
615 {
616 
617 public:
618 
619  //! Empty Constructor
620  DebyeHuckel();
621 
622  //! Copy constructor
623  DebyeHuckel(const DebyeHuckel&);
624 
625  //! Assignment operator
626  DebyeHuckel& operator=(const DebyeHuckel&);
627 
628  //! Full constructor for creating the phase.
629  /*!
630  * @param inputFile File name containing the XML description of the phase
631  * @param id id attribute containing the name of the phase.
632  * (default is the empty string)
633  */
634  DebyeHuckel(std::string inputFile, std::string id = "");
635 
636  //! Full constructor for creating the phase.
637  /*!
638  * @param phaseRef XML phase node containing the description of the phase
639  * @param id id attribute containing the name of the phase.
640  * (default is the empty string)
641  */
642  DebyeHuckel(XML_Node& phaseRef, std::string id = "");
643 
644  /// Destructor.
645  virtual ~DebyeHuckel();
646 
647  //! Duplicator from the ThermoPhase parent class
648  /*!
649  * Given a pointer to a ThermoPhase object, this function will
650  * duplicate the ThermoPhase object and all underlying structures.
651  * This is basically a wrapper around the copy constructor.
652  *
653  * @return returns a pointer to a ThermoPhase
654  */
655  ThermoPhase* duplMyselfAsThermoPhase() const;
656 
657  /**
658  *
659  * @name Utilities
660  * @{
661  */
662 
663  /**
664  * Equation of state type flag. The base class returns
665  * zero. Subclasses should define this to return a unique
666  * non-zero value. Constants defined for this purpose are
667  * listed in mix_defs.h.
668  */
669  virtual int eosType() const;
670 
671  /**
672  * @}
673  * @name Molar Thermodynamic Properties of the Solution --------------
674  * @{
675  */
676 
677  /// Molar enthalpy. Units: J/kmol.
678  /**
679  * Molar enthalpy of the solution. Units: J/kmol.
680  * (HKM -> Bump up to Parent object)
681  */
682  virtual doublereal enthalpy_mole() const;
683 
684  /// Molar internal energy. Units: J/kmol.
685  /**
686  * Molar internal energy of the solution. Units: J/kmol.
687  * (HKM -> Bump up to Parent object)
688  */
689  virtual doublereal intEnergy_mole() const;
690 
691  /// Molar entropy. Units: J/kmol/K.
692  /**
693  * Molar entropy of the solution. Units: J/kmol/K.
694  * For an ideal, constant partial molar volume solution mixture with
695  * pure species phases which exhibit zero volume expansivity:
696  * \f[
697  * \hat s(T, P, X_k) = \sum_k X_k \hat s^0_k(T)
698  * - \hat R \sum_k X_k log(X_k)
699  * \f]
700  * The reference-state pure-species entropies
701  * \f$ \hat s^0_k(T,p_{ref}) \f$ are computed by the
702  * species thermodynamic
703  * property manager. The pure species entropies are independent of
704  * temperature since the volume expansivities are equal to zero.
705  * @see SpeciesThermo
706  *
707  * (HKM -> Bump up to Parent object)
708  */
709  virtual doublereal entropy_mole() const;
710 
711  /// Molar Gibbs function. Units: J/kmol.
712  /*
713  * (HKM -> Bump up to Parent object)
714  */
715  virtual doublereal gibbs_mole() const;
716 
717  /// Molar heat capacity at constant pressure. Units: J/kmol/K.
718  virtual doublereal cp_mole() const;
719 
720  //! Molar heat capacity at constant volume. Units: J/kmol/K.
721  /*
722  * (HKM -> Bump up to Parent object)
723  */
724  virtual doublereal cv_mole() const;
725 
726  //@}
727  /** @name Mechanical Equation of State Properties -------------------------
728  //@{
729  *
730  * In this equation of state implementation, the density is a
731  * function only of the mole fractions. Therefore, it can't be
732  * an independent variable. Instead, the pressure is used as the
733  * independent variable. Functions which try to set the thermodynamic
734  * state by calling setDensity() may cause an exception to be
735  * thrown.
736  */
737 
738  //! Return the thermodynamic pressure (Pa).
739  /*!
740  * For this incompressible system, we return the internally stored
741  * independent value of the pressure.
742  */
743  virtual doublereal pressure() const;
744 
745  //! Set the internally stored pressure (Pa) at constant
746  //! temperature and composition
747  /*!
748  * This method sets the pressure within the object.
749  * The water model is a completely compressible model.
750  * Also, the dielectric constant is pressure dependent.
751  *
752  * @param p input Pressure (Pa)
753  *
754  * @todo Implement a variable pressure capability
755  */
756  virtual void setPressure(doublereal p);
757 
758 protected:
759  //! Calculate the density of the mixture using the partial
760  //! molar volumes and mole fractions as input
761  /*!
762  * The formula for this is
763  *
764  * \f[
765  * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}}
766  * \f]
767  *
768  * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are
769  * the molecular weights, and \f$V_k\f$ are the pure species
770  * molar volumes.
771  *
772  * Note, the basis behind this formula is that in an ideal
773  * solution the partial molar volumes are equal to the pure
774  * species molar volumes. We have additionally specified
775  * in this class that the pure species molar volumes are
776  * independent of temperature and pressure.
777  */
778  virtual void calcDensity();
779 
780 public:
781  //! Set the internally stored density (gm/m^3) of the phase.
782  /*!
783  * Overwritten setDensity() function is necessary because the
784  * density is not an independent variable.
785  *
786  * This function will now throw an error condition
787  *
788  * @internal May have to adjust the strategy here to make
789  * the eos for these materials slightly compressible, in order
790  * to create a condition where the density is a function of
791  * the pressure.
792  *
793  * This function will now throw an error condition if the
794  * input isn't exactly equal to the current density.
795  *
796  *
797  * @todo Now have a compressible ss equation for liquid water.
798  * Therefore, this phase is compressible. May still
799  * want to change the independent variable however.
800  *
801  * NOTE: This is an overwritten function from the State.h
802  * class
803  *
804  * @param rho Input density (kg/m^3).
805  */
806  void setDensity(const doublereal rho);
807 
808  //! Set the internally stored molar density (kmol/m^3) of the phase.
809  /**
810  * Overwritten setMolarDensity() function is necessary because the
811  * density is not an independent variable.
812  *
813  * This function will now throw an error condition if the input
814  * isn't exactly equal to the current molar density.
815  *
816  * NOTE: This is a virtual function overwritten from the State.h
817  * class
818  *
819  * @param conc Input molar density (kmol/m^3).
820  */
821  virtual void setMolarDensity(const doublereal conc);
822 
823  //! Set the temperature (K)
824  /*!
825  * Overwritten setTemperature(double) from State.h. This
826  * function sets the temperature, and makes sure that
827  * the value propagates to underlying objects, such as
828  * the water standard state model.
829  *
830  * @todo Make Phase::setTemperature a virtual function
831  *
832  * @param temp Temperature in kelvin
833  */
834  virtual void setTemperature(const doublereal temp);
835 
836  //! Set the temperature (K) and pressure (Pa)
837  /*!
838  * Set the temperature and pressure.
839  *
840  * @param t Temperature (K)
841  * @param p Pressure (Pa)
842  */
843  virtual void setState_TP(doublereal t, doublereal p);
844 
845  /**
846  * The isothermal compressibility. Units: 1/Pa.
847  * The isothermal compressibility is defined as
848  * \f[
849  * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
850  * \f]
851  */
852  virtual doublereal isothermalCompressibility() const;
853 
854  /**
855  * The thermal expansion coefficient. Units: 1/K.
856  * The thermal expansion coefficient is defined as
857  *
858  * \f[
859  * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
860  * \f]
861  */
862  virtual doublereal thermalExpansionCoeff() const;
863 
864  /**
865  * @}
866  * @name Potential Energy
867  *
868  * Species may have an additional potential energy due to the
869  * presence of external gravitation or electric fields. These
870  * methods allow specifying a potential energy for individual
871  * species.
872  * @{
873  */
874 
875  /**
876  * @}
877  * @name Activities, Standard States, and Activity Concentrations
878  *
879  * The activity \f$a_k\f$ of a species in solution is
880  * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
881  * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
882  * the chemical potential at unit activity, which depends only
883  * on temperature and the pressure.
884  * Activity is assumed to be molality-based here.
885  * @{
886  */
887 
888  //! This method returns an array of generalized concentrations
889  /*!
890  * \f$ C_k\f$ that are defined such that
891  * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
892  * is a standard concentration
893  * defined below. These generalized concentrations are used
894  * by kinetics manager classes to compute the forward and
895  * reverse rates of elementary reactions.
896  *
897  * @param c Array of generalized concentrations. The
898  * units depend upon the implementation of the
899  * reaction rate expressions within the phase.
900  */
901  virtual void getActivityConcentrations(doublereal* c) const;
902 
903  //! Return the standard concentration for the kth species
904  /*!
905  * The standard concentration \f$ C^0_k \f$ used to normalize
906  * the activity (i.e., generalized) concentration in
907  * kinetics calculations.
908  *
909  * For the time being, we will use the concentration of pure
910  * solvent for the the standard concentration of all species.
911  * This has the effect of making reaction rates
912  * based on the molality of species proportional to the
913  * molality of the species.
914  *
915  * @param k Optional parameter indicating the species. The default
916  * is to assume this refers to species 0.
917  * @return
918  * Returns the standard Concentration in units of
919  * m<SUP>3</SUP> kmol<SUP>-1</SUP>.
920  */
921  virtual doublereal standardConcentration(size_t k=0) const;
922 
923  //! Natural logarithm of the standard concentration of the kth species.
924  /*!
925  * @param k index of the species (defaults to zero)
926  */
927  virtual doublereal logStandardConc(size_t k=0) const;
928 
929  //! Returns the units of the standard and generalized concentrations.
930  /*!
931  * Note they have the same units, as their
932  * ratio is defined to be equal to the activity of the kth
933  * species in the solution, which is unitless.
934  *
935  * This routine is used in print out applications where the
936  * units are needed. Usually, MKS units are assumed throughout
937  * the program and in the XML input files.
938  *
939  * The base %ThermoPhase class assigns the default quantities
940  * of (kmol/m3) for all species.
941  * Inherited classes are responsible for overriding the default
942  * values if necessary.
943  *
944  * @param uA Output vector containing the units
945  * uA[0] = kmol units - default = 1
946  * uA[1] = m units - default = -nDim(), the number of spatial
947  * dimensions in the Phase class.
948  * uA[2] = kg units - default = 0;
949  * uA[3] = Pa(pressure) units - default = 0;
950  * uA[4] = Temperature units - default = 0;
951  * uA[5] = time units - default = 0
952  * @param k species index. Defaults to 0.
953  * @param sizeUA output int containing the size of the vector.
954  * Currently, this is equal to 6.
955  */
956  virtual void getUnitsStandardConc(double* uA, int k = 0,
957  int sizeUA = 6) const;
958 
959  //! Get the array of non-dimensional activities at
960  //! the current solution temperature, pressure, and solution concentration.
961  /*!
962  *
963  * We resolve this function at this level by calling
964  * on the activityConcentration function. However,
965  * derived classes may want to override this default
966  * implementation.
967  *
968  * (note solvent is on molar scale).
969  *
970  * @param ac Output vector of activities. Length: m_kk.
971  */
972  virtual void getActivities(doublereal* ac) const;
973 
974  //! Get the array of non-dimensional molality-based
975  //! activity coefficients at
976  //! the current solution temperature, pressure, and solution concentration.
977  /*!
978  * note solvent is on molar scale. The solvent molar
979  * based activity coefficient is returned.
980  *
981  * @param acMolality Vector of Molality-based activity coefficients
982  * Length: m_kk
983  */
984  virtual void
985  getMolalityActivityCoefficients(doublereal* acMolality) const;
986 
987  //@}
988  /// @name Partial Molar Properties of the Solution -----------------
989  //@{
990 
991 
992  //! Get the species chemical potentials. Units: J/kmol.
993  /*!
994  *
995  * This function returns a vector of chemical potentials of the
996  * species in solution.
997  *
998  * \f[
999  * \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} m_k)
1000  * \f]
1001  * or another way to phrase this is
1002  *
1003  * where
1004  *
1005  * @param mu Output vector of species chemical
1006  * potentials. Length: m_kk. Units: J/kmol
1007  */
1008  virtual void getChemPotentials(doublereal* mu) const;
1009 
1010  //! Returns an array of partial molar enthalpies for the species
1011  //! in the mixture. Units (J/kmol)
1012  /*!
1013  * For this phase, the partial molar enthalpies are equal to the
1014  * standard state enthalpies modified by the derivative of the
1015  * molality-based activity coefficient wrt temperature
1016  *
1017  * \f[
1018  * \bar h_k(T,P) = h^{\triangle}_k(T,P) - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
1019  * \f]
1020  * The solvent partial molar enthalpy is equal to
1021  * \f[
1022  * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o}{dT}
1023  * \f]
1024  *
1025  * The temperature dependence of the activity coefficients currently
1026  * only occurs through the temperature dependence of the Debye constant.
1027  *
1028  * @param hbar Output vector of species partial molar enthalpies.
1029  * Length: m_kk. units are J/kmol.
1030  */
1031  virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
1032 
1033  //! Returns an array of partial molar entropies of the species in the
1034  //! solution. Units: J/kmol/K.
1035  /**
1036  * Maxwell's equations provide an insight in how to calculate this
1037  * (p.215 Smith and Van Ness)
1038  *
1039  * d(chemPot_i)/dT = -sbar_i
1040  *
1041  * For this phase, the partial molar entropies are equal to the
1042  * SS species entropies plus the ideal solution contribution.following
1043  * contribution:
1044  * \f[
1045  * \bar s_k(T,P) = \hat s^0_k(T) - R log(M0 * molality[k])
1046  * \f]
1047  * \f[
1048  * \bar s_solvent(T,P) = \hat s^0_solvent(T)
1049  * - R ((xmolSolvent - 1.0) / xmolSolvent)
1050  * \f]
1051  *
1052  * The reference-state pure-species entropies,\f$ \hat s^0_k(T) \f$,
1053  * at the reference pressure, \f$ P_{ref} \f$, are computed by the
1054  * species thermodynamic
1055  * property manager. They are polynomial functions of temperature.
1056  * @see SpeciesThermo
1057  *
1058  * @param sbar Output vector of species partial molar entropies.
1059  * Length = m_kk. units are J/kmol/K.
1060  */
1061  virtual void getPartialMolarEntropies(doublereal* sbar) const;
1062 
1063  //! Return an array of partial molar heat capacities for the
1064  //! species in the mixture. Units: J/kmol/K
1065  /*!
1066  * @param cpbar Output vector of species partial molar heat
1067  * capacities at constant pressure.
1068  * Length = m_kk. units are J/kmol/K.
1069  */
1070  virtual void getPartialMolarCp(doublereal* cpbar) const;
1071 
1072  //! Return an array of partial molar volumes for the
1073  //! species in the mixture. Units: m^3/kmol.
1074  /*!
1075  * For this solution, the partial molar volumes are equal to the
1076  * constant species molar volumes.
1077  *
1078  * @param vbar Output vector of species partial molar volumes.
1079  * Length = m_kk. units are m^3/kmol.
1080  */
1081  virtual void getPartialMolarVolumes(doublereal* vbar) const;
1082 
1083  //@}
1084 
1085 
1086 protected:
1087 
1088  //! Updates the standard state thermodynamic functions at the current T and P of the solution.
1089  /*!
1090  * @internal
1091  *
1092  * This function gets called for every call to a public function in this
1093  * class. It checks to see whether the temperature or pressure has changed and
1094  * thus whether the ss thermodynamics functions must be recalculated.
1095  *
1096  * @param pres Pressure at which to evaluate the standard states.
1097  * The default, indicated by a -1.0, is to use the current pressure
1098  */
1099  //virtual void _updateStandardStateThermo() const;
1100 
1101  //@}
1102  /// @name Thermodynamic Values for the Species Reference States ---
1103  //@{
1104 
1105 
1106  ///////////////////////////////////////////////////////
1107  //
1108  // The methods below are not virtual, and should not
1109  // be overloaded.
1110  //
1111  //////////////////////////////////////////////////////
1112 
1113  /**
1114  * @name Specific Properties
1115  * @{
1116  */
1117 
1118 
1119  /**
1120  * @name Setting the State
1121  *
1122  * These methods set all or part of the thermodynamic
1123  * state.
1124  * @{
1125  */
1126 
1127  //@}
1128 
1129  /**
1130  * @name Chemical Equilibrium
1131  * Chemical equilibrium.
1132  * @{
1133  */
1134 
1135  //!This method is used by the ChemEquil equilibrium solver.
1136  /*!
1137  * It sets the state such that the chemical potentials satisfy
1138  * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
1139  * \left(\frac{\lambda_m} {\hat R T}\right) \f] where
1140  * \f$ \lambda_m \f$ is the element potential of element m. The
1141  * temperature is unchanged. Any phase (ideal or not) that
1142  * implements this method can be equilibrated by ChemEquil.
1143  *
1144  * @param lambda_RT Input vector of dimensionless element potentials
1145  * The length is equal to nElements().
1146  */
1147 public:
1148  virtual void setToEquilState(const doublereal* lambda_RT) {
1149  err("setToEquilState");
1150  }
1151 
1152 
1153  //@}
1154 
1155 
1156  //! Set the equation of state parameters
1157  /*!
1158  * @internal
1159  * The number and meaning of these depends on the subclass.
1160  *
1161  * @param n number of parameters
1162  * @param c array of \a n coefficients
1163  */
1164  virtual void setParameters(int n, doublereal* const c);
1165 
1166  //! Get the equation of state parameters in a vector
1167  /*!
1168  * @internal
1169  * The number and meaning of these depends on the subclass.
1170  *
1171  * @param n number of parameters
1172  * @param c array of \a n coefficients
1173  */
1174  virtual void getParameters(int& n, doublereal* const c) const;
1175 
1176  //! Set equation of state parameter values from XML entries.
1177  /*!
1178  *
1179  * This method is called by function importPhase() in
1180  * file importCTML.cpp when processing a phase definition in
1181  * an input file. It should be overloaded in subclasses to set
1182  * any parameters that are specific to that particular phase
1183  * model. Note, this method is called before the phase is
1184  * initialized with elements and/or species.
1185  *
1186  * @param eosdata An XML_Node object corresponding to
1187  * the "thermo" entry for this phase in the input file.
1188  */
1189  virtual void setParametersFromXML(const XML_Node& eosdata);
1190 
1191  //---------------------------------------------------------
1192  /// @name Critical state properties.
1193  /// These methods are only implemented by some subclasses.
1194 
1195  //@{
1196 
1197 
1198 
1199  //@}
1200 
1201  /// @name Saturation properties.
1202  /// These methods are only implemented by subclasses that
1203  /// implement full liquid-vapor equations of state.
1204  ///
1205  virtual doublereal satTemperature(doublereal p) const {
1206  err("satTemperature");
1207  return -1.0;
1208  }
1209 
1210  //! Get the saturation pressure for a given temperature.
1211  /*!
1212  * Note the limitations of this function. Stability considerations
1213  * concerning multiphase equilibrium are ignored in this
1214  * calculation. Therefore, the call is made directly to the SS of
1215  * water underneath. The object is put back into its original
1216  * state at the end of the call.
1217  *
1218  * @todo This is probably not implemented correctly. The stability
1219  * of the salt should be added into this calculation. The
1220  * underlying water model may be called to get the stability
1221  * of the pure water solution, if needed.
1222  *
1223  * @param T Temperature (kelvin)
1224  */
1225  virtual doublereal satPressure(doublereal T) const {
1226  err("satPressure");
1227  return -1.0;
1228  }
1229 
1230  virtual doublereal vaporFraction() const {
1231  err("vaprFraction");
1232  return -1.0;
1233  }
1234 
1235  virtual void setState_Tsat(doublereal t, doublereal x) {
1236  err("setState_sat");
1237  }
1238 
1239  virtual void setState_Psat(doublereal p, doublereal x) {
1240  err("setState_sat");
1241  }
1242 
1243  //@}
1244 
1245 
1246  /*
1247  * -------------- Utilities -------------------------------
1248  */
1249 
1250 
1251  //! Initialize the object's internal lengths after species are set
1252  /**
1253  * @internal Initialize. This method is provided to allow
1254  * subclasses to perform any initialization required after all
1255  * species have been added. For example, it might be used to
1256  * resize internal work arrays that must have an entry for
1257  * each species. The base class implementation does nothing,
1258  * and subclasses that do not require initialization do not
1259  * need to overload this method. When importing a CTML phase
1260  * description, this method is called just prior to returning
1261  * from function importPhase().
1262  *
1263  * Cascading call sequence downwards starting with Parent.
1264  *
1265  * @internal
1266  *
1267  * @see importCTML.cpp
1268  */
1269  virtual void initThermo();
1270 
1271 
1272  //! Initialization of a DebyeHuckel phase using an xml file
1273  /*!
1274  * This routine is a precursor to initThermo(XML_Node*)
1275  * routine, which does most of the work.
1276  *
1277  * @param infile XML file containing the description of the
1278  * phase
1279  *
1280  * @param id Optional parameter identifying the name of the
1281  * phase. If none is given, the first XML
1282  * phase element will be used.
1283  */
1284  void constructPhaseFile(std::string infile, std::string id="");
1285 
1286  //! Import and initialize a DebyeHuckel phase
1287  //! specification in an XML tree into the current object.
1288  /*!
1289  * Here we read an XML description of the phase.
1290  * We import descriptions of the elements that make up the
1291  * species in a phase.
1292  * We import information about the species, including their
1293  * reference state thermodynamic polynomials. We then freeze
1294  * the state of the species.
1295  *
1296  * Then, we read the species molar volumes from the xml
1297  * tree to finish the initialization.
1298  *
1299  * @param phaseNode This object must be the phase node of a
1300  * complete XML tree
1301  * description of the phase, including all of the
1302  * species data. In other words while "phase" must
1303  * point to an XML phase object, it must have
1304  * sibling nodes "speciesData" that describe
1305  * the species in the phase.
1306  *
1307  * @param id ID of the phase. If nonnull, a check is done
1308  * to see if phaseNode is pointing to the phase
1309  * with the correct id.
1310  */
1311  void constructPhaseXML(XML_Node& phaseNode, std::string id="");
1312 
1313  //! Process the XML file after species are set up.
1314  /*!
1315  * This gets called from importPhase(). It processes the XML file
1316  * after the species are set up. This is the main routine for
1317  * reading in activity coefficient parameters.
1318  *
1319  * @param phaseNode This object must be the phase node of a
1320  * complete XML tree
1321  * description of the phase, including all of the
1322  * species data. In other words while "phase" must
1323  * point to an XML phase object, it must have
1324  * sibling nodes "speciesData" that describe
1325  * the species in the phase.
1326  * @param id ID of the phase. If nonnull, a check is done
1327  * to see if phaseNode is pointing to the phase
1328  * with the correct id.
1329  */
1330  virtual void initThermoXML(XML_Node& phaseNode, std::string id);
1331 
1332  //! Return the Debye Huckel constant as a function of temperature
1333  //! and pressure (Units = sqrt(kg/gmol))
1334  /*!
1335  * The default is to assume that it is constant, given
1336  * in the initialization process, and stored in the
1337  * member double, m_A_Debye. Optionally, a full water treatment may be employed that makes
1338  * \f$ A_{Debye} \f$ a full function of T and P.
1339  *
1340  * \f[
1341  * A_{Debye} = \frac{F e B_{Debye}}{8 \pi \epsilon R T} {\left( C_o \tilde{M}_o \right)}^{1/2}
1342  * \f]
1343  * where
1344  *
1345  * \f[
1346  * B_{Debye} = \frac{F} {{(\frac{\epsilon R T}{2})}^{1/2}}
1347  * \f]
1348  * Therefore:
1349  * \f[
1350  * A_{Debye} = \frac{1}{8 \pi}
1351  * {\left(\frac{2 N_a \rho_o}{1000}\right)}^{1/2}
1352  * {\left(\frac{N_a e^2}{\epsilon R T }\right)}^{3/2}
1353  * \f]
1354  *
1355  * Units = sqrt(kg/gmol)
1356  *
1357  * where
1358  * - \f$ N_a \f$ is Avogadro's number
1359  * - \f$ \rho_w \f$ is the density of water
1360  * - \f$ e \f$ is the electronic charge
1361  * - \f$ \epsilon = K \epsilon_o \f$ is the permittivity of water
1362  * where \f$ K \f$ is the dielectric constant of water,
1363  * and \f$ \epsilon_o \f$ is the permittivity of free space.
1364  * = \f$ \rho_o \f$ is the density of the solvent in its standard state.
1365  *
1366  * Nominal value at 298 K and 1 atm = 1.172576 (kg/gmol)<SUP>1/2</SUP>
1367  * based on:
1368  * - \f$ \epsilon / \epsilon_0 \f$ = 78.54
1369  * (water at 25C)
1370  * - T = 298.15 K
1371  * - B_Debye = 3.28640E9 (kg/gmol)<SUP>1/2</SUP> m<SUP>-1</SUP>
1372  *
1373  * @param temperature Temperature in kelvin. Defaults to -1, in which
1374  * case the temperature of the phase is assumed.
1375  *
1376  * @param pressure Pressure (Pa). Defaults to -1, in which
1377  * case the pressure of the phase is assumed.
1378  */
1379  virtual double A_Debye_TP(double temperature = -1.0,
1380  double pressure = -1.0) const;
1381 
1382 
1383  //! Value of the derivative of the Debye Huckel constant with
1384  //! respect to temperature.
1385  /*!
1386  * This is a function of temperature and pressure. See A_Debye_TP() for
1387  * a definition of \f$ A_{Debye} \f$.
1388  *
1389  * Units = sqrt(kg/gmol) K-1
1390  *
1391  * @param temperature Temperature in kelvin. Defaults to -1, in which
1392  * case the temperature of the phase is assumed.
1393  *
1394  * @param pressure Pressure (Pa). Defaults to -1, in which
1395  * case the pressure of the phase is assumed.
1396  */
1397  virtual double dA_DebyedT_TP(double temperature = -1.0,
1398  double pressure = -1.0) const;
1399 
1400  //! Value of the 2nd derivative of the Debye Huckel constant with
1401  //! respect to temperature as a function of temperature and pressure.
1402  /*!
1403  * This is a function of temperature and pressure. See A_Debye_TP() for
1404  * a definition of \f$ A_{Debye} \f$.
1405  *
1406  * Units = sqrt(kg/gmol) K-2
1407  *
1408  * @param temperature Temperature in kelvin. Defaults to -1, in which
1409  * case the temperature of the phase is assumed.
1410  *
1411  * @param pressure Pressure (Pa). Defaults to -1, in which
1412  * case the pressure of the phase is assumed.
1413  */
1414  virtual double d2A_DebyedT2_TP(double temperature = -1.0,
1415  double pressure = -1.0) const;
1416 
1417  //! Value of the derivative of the Debye Huckel constant with
1418  //! respect to pressure, as a function of temperature and pressure.
1419  /*!
1420  * This is a function of temperature and pressure. See A_Debye_TP() for
1421  * a definition of \f$ A_{Debye} \f$.
1422  *
1423  * Units = sqrt(kg/gmol) Pa-1
1424  *
1425  * @param temperature Temperature in kelvin. Defaults to -1, in which
1426  * case the temperature of the phase is assumed.
1427  *
1428  * @param pressure Pressure (Pa). Defaults to -1, in which
1429  * case the pressure of the phase is assumed.
1430  */
1431  virtual double dA_DebyedP_TP(double temperature = -1.0,
1432  double pressure = -1.0) const;
1433 
1434  //!Reports the ionic radius of the kth species
1435  /*!
1436  * @param k species index.
1437  */
1438  double AionicRadius(int k = 0) const;
1439 
1440  //! Returns the form of the Debye-Huckel parameterization used
1441  int formDH() const {
1442  return m_formDH;
1443  }
1444 
1445  //! Returns a reference to M_Beta_ij
1446  Array2D& get_Beta_ij() {
1447  return m_Beta_ij;
1448  }
1449 
1450 private:
1451 
1452 
1453  //! Static function that implements the non-polar species
1454  //! salt-out modifications.
1455  /*!
1456  * Returns the calculated activity coefficients.
1457  *
1458  * @param IionicMolality Value of the ionic molality (sqrt(gmol/kg))
1459  */
1460  double _nonpolarActCoeff(double IionicMolality) const;
1461 
1462 
1463  //! Formula for the osmotic coefficient that occurs in the GWB.
1464  /*!
1465  * It is originally from Helgeson for a variable
1466  * NaCl brine. It's to be used with extreme caution.
1467  */
1468  double _osmoticCoeffHelgesonFixedForm() const;
1469 
1470  //! Formula for the log of the water activity that occurs in the GWB.
1471  /*!
1472  * It is originally from Helgeson for a variable
1473  * NaCl brine. It's to be used with extreme caution.
1474  */
1475  double _lnactivityWaterHelgesonFixedForm() const;
1476 
1477 
1478  //@}
1479 
1480 
1481 protected:
1482 
1483 
1484  //! form of the Debye-Huckel parameterization used in the model.
1485  /*!
1486  * The options are described at the top of this document,
1487  * and in the general documentation.
1488  * The list is repeated here:
1489  *
1490  * DHFORM_DILUTE_LIMIT = 0 (default)
1491  * DHFORM_BDOT_AK = 1
1492  * DHFORM_BDOT_AUNIFORM = 2
1493  * DHFORM_BETAIJ = 3
1494  * DHFORM_PITZER_BETAIJ = 4
1495  */
1496  int m_formDH;
1497 
1498  /**
1499  * Format for the generalized concentration:
1500  *
1501  * 0 = unity
1502  * 1 = molar_volume
1503  * 2 = solvent_volume (default)
1504  *
1505  * The generalized concentrations can have three different forms
1506  * depending on the value of the member attribute m_formGC, which
1507  * is supplied in the constructor.
1508  * <TABLE>
1509  * <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
1510  * <TR><TD> 0 </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
1511  * <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
1512  * <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
1513  * </TABLE>
1514  *
1515  * The value and form of the generalized concentration will affect
1516  * reaction rate constants involving species in this phase.
1517  *
1518  * (HKM Note: Using option #1 may lead to spurious results and
1519  * has been included only with warnings. The reason is that it
1520  * molar volumes of electrolytes may often be negative. The
1521  * molar volume of H+ is defined to be zero too. Either options
1522  * 0 or 2 are the appropriate choice. Option 0 leads to
1523  * bulk reaction rate constants which have units of s-1.
1524  * Option 2 leads to bulk reaction rate constants for
1525  * bimolecular rxns which have units of m-3 kmol-1 s-1.)
1526  */
1527  int m_formGC;
1528 
1529  //! Vector containing the electrolyte species type
1530  /*!
1531  * The possible types are:
1532  * - solvent
1533  * - Charged Species
1534  * - weakAcidAssociated
1535  * - strongAcidAssociated
1536  * - polarNeutral
1537  * - nonpolarNeutral
1538  * .
1539  */
1540  vector_int m_electrolyteSpeciesType;
1541 
1542  /**
1543  * a_k = Size of the ionic species in the DH formulation
1544  * units = meters
1545  */
1546  vector_fp m_Aionic;
1547 
1548  /**
1549  * Current value of the ionic strength on the molality scale
1550  */
1551  mutable double m_IionicMolality;
1552 
1553  /**
1554  * Maximum value of the ionic strength allowed in the
1555  * calculation of the activity coefficients.
1556  */
1557  double m_maxIionicStrength;
1558 
1559 public:
1560 
1561  /**
1562  * If true, then the fixed for of Helgeson's activity
1563  * for water is used instead of the rigorous form
1564  * obtained from Gibbs-Duhem relation. This should be
1565  * used with caution, and is really only included as a
1566  * validation exercise.
1567  */
1568  bool m_useHelgesonFixedForm;
1569 protected:
1570 
1571  //! Stoichiometric ionic strength on the molality scale
1572  mutable double m_IionicMolalityStoich;
1573 
1574 public:
1575 
1576  /**
1577  * Form of the constant outside the Debye-Huckel term
1578  * called A. It's normally a function of temperature
1579  * and pressure. However, it can be set from the
1580  * input file in order to aid in numerical comparisons.
1581  * Acceptable forms:
1582  *
1583  * A_DEBYE_CONST 0
1584  * A_DEBYE_WATER 1
1585  *
1586  * The A_DEBYE_WATER form may be used for water solvents
1587  * with needs to cover varying temperatures and pressures.
1588  * Note, the dielectric constant of water is a relatively
1589  * strong function of T, and its variability must be
1590  * accounted for,
1591  */
1592  int m_form_A_Debye;
1593 
1594 protected:
1595 
1596  //! Current value of the Debye Constant, A_Debye
1597  /**
1598  * A_Debye -> this expression appears on the top of the
1599  * ln actCoeff term in the general Debye-Huckel
1600  * expression
1601  * It depends on temperature and pressure.
1602  *
1603  * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
1604  *
1605  * Units = sqrt(kg/gmol)
1606  *
1607  * Nominal value(298K, atm) = 1.172576 sqrt(kg/gmol)
1608  * based on:
1609  * epsilon/epsilon_0 = 78.54
1610  * (water at 25C)
1611  * T = 298.15 K
1612  * B_Debye = 3.28640E9 sqrt(kg/gmol)/m
1613  *
1614  * note in Pitzer's nomenclature, A_phi = A_Debye/3.0
1615  */
1616  mutable double m_A_Debye;
1617 
1618  //! Current value of the constant that appears in the denominator
1619  /**
1620  * B_Debye -> this expression appears on the bottom of the
1621  * ln actCoeff term in the general Debye-Huckel
1622  * expression
1623  * It depends on temperature
1624  *
1625  * B_Bebye = F / sqrt( epsilon R T / 2 )
1626  *
1627  * Units = sqrt(kg/gmol) / m
1628  *
1629  * Nominal value = 3.28640E9 sqrt(kg/gmol) / m
1630  * based on:
1631  * epsilon/epsilon_0 = 78.54
1632  * (water at 25C)
1633  * T = 298.15 K
1634  */
1635  double m_B_Debye;
1636 
1637  //! Array of B_Dot values
1638  /**
1639  * B_Dot -> This expression is an extension of the
1640  * Debye-Huckel expression used in some formulations
1641  * to extend DH to higher molalities.
1642  * B_dot is specific to the major ionic pair.
1643  */
1644  vector_fp m_B_Dot;
1645 
1646  /**
1647  * m_npActCoeff -> These are coefficients to describe
1648  * the increase in activity coeff for non-polar molecules
1649  * due to the electrolyte becoming stronger (the so-called
1650  * salt-out effect)
1651  */
1652  vector_fp m_npActCoeff;
1653 
1654 
1655  //! Pointer to the Water standard state object
1656  /*!
1657  * derived from the equation of state for water.
1658  */
1659  PDSS_Water* m_waterSS;
1660 
1661  //! Storage for the density of water's standard state
1662  /*!
1663  * Density depends on temperature and pressure.
1664  */
1665  double m_densWaterSS;
1666 
1667  /**
1668  * Pointer to the water property calculator
1669  */
1670  WaterProps* m_waterProps;
1671 
1672  /**
1673  * Vector containing the species reference exp(-G/RT) functions
1674  * at T = m_tlast
1675  */
1676  mutable vector_fp m_expg0_RT;
1677 
1678  /**
1679  * Vector of potential energies for the species.
1680  */
1681  mutable vector_fp m_pe;
1682 
1683  /**
1684  * Temporary array used in equilibrium calculations
1685  */
1686  mutable vector_fp m_pp;
1687 
1688  /**
1689  * vector of size m_kk, used as a temporary holding area.
1690  */
1691  mutable vector_fp m_tmpV;
1692 
1693  /**
1694  * Stoichiometric species charge -> This is for calculations
1695  * of the ionic strength which ignore ion-ion pairing into
1696  * neutral molecules. The Stoichiometric species charge is the
1697  * charge of one of the ion that would occur if the species broke
1698  * into two charged ion pairs.
1699  * NaCl -> m_speciesCharge_Stoich = -1;
1700  * HSO4- -> H+ + SO42- = -2
1701  * -> The other charge is calculated.
1702  * For species that aren't ion pairs, it's equal to the
1703  * m_speciesCharge[] value.
1704  */
1705  vector_fp m_speciesCharge_Stoich;
1706 
1707  /**
1708  * Array of 2D data used in the DHFORM_BETAIJ formulation
1709  * Beta_ij.value(i,j) is the coefficient of the jth species
1710  * for the specification of the chemical potential of the ith
1711  * species.
1712  */
1713  Array2D m_Beta_ij;
1714 
1715  //! Logarithm of the activity coefficients on the molality scale.
1716  /*!
1717  * mutable because we change this if the composition
1718  * or temperature or pressure changes.
1719  */
1720  mutable vector_fp m_lnActCoeffMolal;
1721 
1722  //! Derivative of log act coeff wrt T
1723  mutable vector_fp m_dlnActCoeffMolaldT;
1724 
1725  //! 2nd Derivative of log act coeff wrt T
1726  mutable vector_fp m_d2lnActCoeffMolaldT2;
1727 
1728  //! Derivative of log act coeff wrt P
1729  mutable vector_fp m_dlnActCoeffMolaldP;
1730 
1731 private:
1732  doublereal err(std::string msg) const;
1733 
1734  //! Initialize the internal lengths.
1735  /*!
1736  * This internal function adjusts the lengths of arrays based on
1737  * the number of species.
1738  */
1739  void initLengths();
1740 
1741 private:
1742  //! Calculate the log activity coefficients
1743  /*!
1744  * This function updates the internally stored
1745  * natural logarithm of the molality activity coefficients
1746  */
1747  void s_update_lnMolalityActCoeff() const;
1748 
1749  //! Calculation of temperature derivative of activity coefficient
1750  /*!
1751  * Using internally stored values, this function calculates
1752  * the temperature derivative of the logarithm of the
1753  * activity coefficient for all species in the mechanism.
1754  *
1755  * We assume that the activity coefficients are current in this routine
1756  *
1757  *
1758  *
1759  * The solvent activity coefficient is on the molality scale. Its derivative is too.
1760  */
1761  void s_update_dlnMolalityActCoeff_dT() const;
1762 
1763  //! Calculate the temperature 2nd derivative of the activity coefficient
1764  /*!
1765  * Using internally stored values, this function calculates
1766  * the temperature 2nd derivative of the logarithm of the
1767  * activity coefficient for all species in the mechanism.
1768  *
1769  * We assume that the activity coefficients are current in this routine
1770  *
1771  * solvent activity coefficient is on the molality
1772  * scale. Its derivatives are too.
1773  *
1774  * note: private routine
1775  */
1776  void s_update_d2lnMolalityActCoeff_dT2() const;
1777 
1778  //! Calculate the pressure derivative of the activity coefficient
1779  /*!
1780  * Using internally stored values, this function calculates
1781  * the pressure derivative of the logarithm of the
1782  * activity coefficient for all species in the mechanism.
1783  *
1784  * We assume that the activity coefficients, molalities,
1785  * and A_Debye are current.
1786  *
1787  * solvent activity coefficient is on the molality
1788  * scale. Its derivatives are too.
1789  */
1790  void s_update_dlnMolalityActCoeff_dP() const;
1791 };
1792 
1793 }
1794 
1795 #endif
1796 
1797 
1798 
1799 
1800 
Generated by   doxygen 1.8.2