Loading [MathJax]/extensions/tex2jax.js
Cantera  3.2.0a1
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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// This file is part of Cantera. See License.txt in the top-level directory or
12// at https://cantera.org/license.txt for license and copyright information.
13
14#ifndef CT_DEBYEHUCKEL_H
15#define CT_DEBYEHUCKEL_H
16
17#include "MolalityVPSSTP.h"
18#include "cantera/base/Array.h"
19
20namespace Cantera
21{
22
23//! @name Formats for the Activity Coefficients
24//!
25//! These are possible formats for the molality-based activity coefficients.
26//! @{
27
28#define DHFORM_DILUTE_LIMIT 0
29#define DHFORM_BDOT_AK 1
30#define DHFORM_BDOT_ACOMMON 2
31#define DHFORM_BETAIJ 3
32#define DHFORM_PITZER_BETAIJ 4
33
34//! @}
35//! @name Acceptable ways to calculate the value of A_Debye
36//! @{
37
38#define A_DEBYE_CONST 0
39#define A_DEBYE_WATER 1
40//! @}
41
42class WaterProps;
43class PDSS_Water;
44
45/**
46 * @ingroup thermoprops
47 *
48 * Class DebyeHuckel represents a dilute liquid electrolyte phase which obeys
49 * the Debye Huckel formulation for nonideality.
50 *
51 * The concentrations of the ionic species are assumed to obey the
52 * electroneutrality condition.
53 *
54 * ## Specification of Species Standard State Properties
55 *
56 * The standard states are on the unit molality basis. Therefore, in the
57 * documentation below, the normal @f$ o @f$ superscript is replaced with the
58 * @f$ \triangle @f$ symbol. The reference state symbol is now
59 * @f$ \triangle, ref @f$.
60 *
61 * It is assumed that the reference state thermodynamics may be obtained by a
62 * pointer to a populated species thermodynamic property manager class (see
63 * ThermoPhase::m_spthermo). How to relate pressure changes to the reference
64 * state thermodynamics is resolved at this level.
65 *
66 * For an incompressible, stoichiometric substance, the molar internal energy is
67 * independent of pressure. Since the thermodynamic properties are specified by
68 * giving the standard-state enthalpy, the term @f$ P_0 \hat v @f$ is subtracted
69 * from the specified molar enthalpy to compute the molar internal energy. The
70 * entropy is assumed to be independent of the pressure.
71 *
72 * The enthalpy function is given by the following relation.
73 *
74 * @f[
75 * h^\triangle_k(T,P) = h^{\triangle,ref}_k(T)
76 * + \tilde v \left( P - P_{ref} \right)
77 * @f]
78 *
79 * For an incompressible, stoichiometric substance, the molar internal energy is
80 * independent of pressure. Since the thermodynamic properties are specified by
81 * giving the standard-state enthalpy, the term @f$ P_{ref} \tilde v @f$ is
82 * subtracted from the specified reference molar enthalpy to compute the molar
83 * internal energy.
84 *
85 * @f[
86 * u^\triangle_k(T,P) = h^{\triangle,ref}_k(T) - P_{ref} \tilde v
87 * @f]
88 *
89 * The standard state heat capacity and entropy are independent of pressure. The
90 * standard state Gibbs free energy is obtained from the enthalpy and entropy
91 * functions.
92 *
93 * The current model assumes that an incompressible molar volume for all
94 * solutes. The molar volume for the water solvent, however, is obtained from a
95 * pure water equation of state, waterSS. Therefore, the water standard state
96 * varies with both T and P. It is an error to request standard state water
97 * properties at a T and P where the water phase is not a stable phase, that is,
98 * beyond its spinodal curve.
99 *
100 * ## Specification of Solution Thermodynamic Properties
101 *
102 * Chemical potentials of the solutes, @f$ \mu_k @f$, and the solvent, @f$ \mu_o
103 * @f$, which are based on the molality form, have the following general format:
104 *
105 * @f[
106 * \mu_k = \mu^{\triangle}_k(T,P) + R T \ln(\gamma_k^{\triangle} \frac{m_k}{m^\triangle})
107 * @f]
108 * @f[
109 * \mu_o = \mu^o_o(T,P) + RT \ln(a_o)
110 * @f]
111 *
112 * where @f$ \gamma_k^{\triangle} @f$ is the molality based activity coefficient
113 * for species @f$ k @f$.
114 *
115 * Individual activity coefficients of ions can not be independently measured.
116 * Instead, only binary pairs forming electroneutral solutions can be measured.
117 *
118 * ### Ionic Strength
119 *
120 * Most of the parameterizations within the model use the ionic strength as a
121 * key variable. The ionic strength, @f$ I @f$ is defined as follows
122 *
123 * @f[
124 * I = \frac{1}{2} \sum_k{m_k z_k^2}
125 * @f]
126 *
127 * @f$ m_k @f$ is the molality of the kth species. @f$ z_k @f$ is the charge of
128 * the kth species. Note, the ionic strength is a defined units quantity. The
129 * molality has defined units of gmol kg-1, and therefore the ionic strength has
130 * units of sqrt( gmol kg-1).
131 *
132 * In some instances, from some authors, a different formulation is used for the
133 * ionic strength in the equations below. The different formulation is due to
134 * the possibility of the existence of weak acids and how association wrt to the
135 * weak acid equilibrium relation affects the calculation of the activity
136 * coefficients via the assumed value of the ionic strength.
137 *
138 * If we are to assume that the association reaction doesn't have an effect on
139 * the ionic strength, then we will want to consider the associated weak acid as
140 * in effect being fully dissociated, when we calculate an effective value for
141 * the ionic strength. We will call this calculated value, the stoichiometric
142 * ionic strength, @f$ I_s @f$, putting a subscript s to denote it from the more
143 * straightforward calculation of @f$ I @f$.
144 *
145 * @f[
146 * I_s = \frac{1}{2} \sum_k{m_k^s z_k^2}
147 * @f]
148 *
149 * Here, @f$ m_k^s @f$ is the value of the molalities calculated assuming that
150 * all weak acid-base pairs are in their fully dissociated states. This
151 * calculation may be simplified by considering that the weakly associated acid
152 * may be made up of two charged species, k1 and k2, each with their own
153 * charges, obeying the following relationship:
154 *
155 * @f[
156 * z_k = z_{k1} + z_{k2}
157 * @f]
158 * Then, we may only need to specify one charge value, say, @f$ z_{k1} @f$, the
159 * cation charge number, in order to get both numbers, since we have already
160 * specified @f$ z_k @f$ in the definition of original species. Then, the
161 * stoichiometric ionic strength may be calculated via the following formula.
162 *
163 * @f[
164 * I_s = \frac{1}{2} \left(\sum_{k,ions}{m_k z_k^2}+
165 * \sum_{k,weak_assoc}(m_k z_{k1}^2 + m_k z_{k2}^2) \right)
166 * @f]
167 *
168 * The specification of which species are weakly associated acids is made in YAML
169 * input files by specifying the corresponding charge @f$ k1 @f$ as the `weak-acid-charge`
170 * parameter of the `Debye-Huckel` block in the corresponding species entry.
171 *
172 * Because we need the concept of a weakly associated acid in order to calculate
173 * @f$ I_s @f$ we need to catalog all species in the phase. This is done using
174 * the following categories:
175 *
176 * - `cEST_solvent` Solvent species (neutral)
177 * - `cEST_chargedSpecies` Charged species (charged)
178 * - `cEST_weakAcidAssociated` Species which can break apart into charged species.
179 * It may or may not be charged. These may or
180 * may not be be included in the
181 * species solution vector.
182 * - `cEST_strongAcidAssociated` Species which always breaks apart into charged species.
183 * It may or may not be charged. Normally, these aren't included
184 * in the speciation vector.
185 * - `cEST_polarNeutral` Polar neutral species
186 * - `cEST_nonpolarNeutral` Non polar neutral species
187 *
188 * Polar and non-polar neutral species are differentiated, because some
189 * additions to the activity coefficient expressions distinguish between these
190 * two types of solutes. This is the so-called salt-out effect.
191 *
192 * In a YAML input file, the type of species is specified in the
193 * `electrolyte-species-type` field of the `Debye-Huckel` block in the corresponding
194 * species entry. Note, this is not considered a part of the specification of the
195 * standard state for the species, at this time.
196 *
197 * Much of the species electrolyte type information is inferred from other
198 * information in the input file. For example, as species which is charged is
199 * given the "chargedSpecies" default category. A neutral solute species is put
200 * into the "nonpolarNeutral" category by default.
201 *
202 * The specification of solute activity coefficients depends on the model
203 * assumed for the Debye-Huckel term. The model is set by the internal parameter
204 * #m_formDH. We will now describe each category in its own section.
205 *
206 * ### Debye-Huckel Dilute Limit
207 *
208 * DHFORM_DILUTE_LIMIT = 0
209 *
210 * This form assumes a dilute limit to DH, and is mainly for informational purposes:
211 * @f[
212 * \ln(\gamma_k^\triangle) = - z_k^2 A_{Debye} \sqrt{I}
213 * @f]
214 * where @f$ I @f$ is the ionic strength
215 * @f[
216 * I = \frac{1}{2} \sum_k{m_k z_k^2}
217 * @f]
218 *
219 * The activity for the solvent water,@f$ a_o @f$, is not independent and must
220 * be determined from the Gibbs-Duhem relation.
221 *
222 * @f[
223 * \ln(a_o) = \frac{X_o - 1.0}{X_o} + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{3/2}
224 * @f]
225 *
226 * ### Bdot Formulation
227 *
228 * DHFORM_BDOT_AK = 1
229 *
230 * This form assumes Bethke's format for the Debye Huckel activity coefficient:
231 *
232 * @f[
233 * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye} \sqrt{I}}{ 1 + B_{Debye} a_k \sqrt{I}}
234 * + \ln(10) B^{dot}_k I
235 * @f]
236 *
237 * Note, this particular form where @f$ a_k @f$ can differ in multielectrolyte
238 * solutions has problems with respect to a Gibbs-Duhem analysis. However, we
239 * include it here because there is a lot of data fit to it.
240 *
241 * The activity for the solvent water,@f$ a_o @f$, is not independent and must
242 * be determined from the Gibbs-Duhem relation. Here, we use:
243 *
244 * @f[
245 * \ln(a_o) = \frac{X_o - 1.0}{X_o}
246 * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{1/2}
247 * \left[ \sum_k{\frac{1}{2} m_k z_k^2 \sigma( B_{Debye} a_k \sqrt{I} ) } \right]
248 * - \frac{\ln 10}{2} \tilde{M}_o I \sum_k{ B^{dot}_k m_k}
249 * @f]
250 * where
251 * @f[
252 * \sigma (y) = \frac{3}{y^3} \left[ (1+y) - 2 \ln(1 + y) - \frac{1}{1+y} \right]
253 * @f]
254 *
255 * Additionally, Helgeson's formulation for the water activity is offered as an
256 * alternative.
257 *
258 * ### Bdot Formulation with Uniform Size Parameter in the Denominator
259 *
260 * DHFORM_BDOT_AUNIFORM = 2
261 *
262 * This form assumes Bethke's format for the Debye-Huckel activity coefficient
263 *
264 * @f[
265 * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye} \sqrt{I}}{ 1 + B_{Debye} a \sqrt{I}}
266 * + \ln(10) B^{dot}_k I
267 * @f]
268 *
269 * The value of a is determined at the beginning of the calculation, and not changed.
270 *
271 * @f[
272 * \ln(a_o) = \frac{X_o - 1.0}{X_o}
273 * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{3/2} \sigma( B_{Debye} a \sqrt{I} )
274 * - \frac{\ln 10}{2} \tilde{M}_o I \sum_k{ B^{dot}_k m_k}
275 * @f]
276 *
277 * ### Beta_IJ formulation
278 *
279 * DHFORM_BETAIJ = 3
280 *
281 * This form assumes a linear expansion in a virial coefficient form. It is used
282 * extensively in the book by Newmann, "Electrochemistry Systems", and is the
283 * beginning of more complex treatments for stronger electrolytes, fom Pitzer
284 * and from Harvey, Moller, and Weire.
285 *
286 * @f[
287 * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye} \sqrt{I}}{ 1 + B_{Debye} a \sqrt{I}}
288 * + 2 \sum_j \beta_{j,k} m_j
289 * @f]
290 *
291 * In the current treatment the binary interaction coefficients, @f$
292 * \beta_{j,k} @f$, are independent of temperature and pressure.
293 *
294 * @f[
295 * \ln(a_o) = \frac{X_o - 1.0}{X_o}
296 * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} (I)^{3/2} \sigma( B_{Debye} a \sqrt{I} )
297 * - \tilde{M}_o \sum_j \sum_k \beta_{j,k} m_j m_k
298 * @f]
299 *
300 * In this formulation the ionic radius, @f$ a @f$, is a constant, specified as part
301 * of the species definition.
302 *
303 * The @f$ \beta_{j,k} @f$ parameters are binary interaction parameters. There are in
304 * principle @f$ N (N-1) /2 @f$ different, symmetric interaction parameters,
305 * where @f$ N @f$ are the number of solute species in the mechanism.
306 *
307 * ### Pitzer Beta_IJ formulation
308 *
309 * DHFORM_PITZER_BETAIJ = 4
310 *
311 * This form assumes an activity coefficient formulation consistent with a
312 * truncated form of Pitzer's formulation. Pitzer's formulation is equivalent to
313 * the formulations above in the dilute limit, where rigorous theory may be
314 * applied.
315 *
316 * @f[
317 * \ln(\gamma_k^\triangle) = -z_k^2 \frac{A_{Debye}}{3} \frac{\sqrt{I}}{ 1 + B_{Debye} a \sqrt{I}}
318 * -2 z_k^2 \frac{A_{Debye}}{3} \frac{\ln(1 + B_{Debye} a \sqrt{I})}{ B_{Debye} a}
319 * + 2 \sum_j \beta_{j,k} m_j
320 * @f]
321 * @f[
322 * \ln(a_o) = \frac{X_o - 1.0}{X_o}
323 * + \frac{ 2 A_{Debye} \tilde{M}_o}{3} \frac{(I)^{3/2} }{1 + B_{Debye} a \sqrt{I} }
324 * - \tilde{M}_o \sum_j \sum_k \beta_{j,k} m_j m_k
325 * @f]
326 *
327 * ### Specification of the Debye Huckel Constants
328 *
329 * In the equations above, the formulas for @f$ A_{Debye} @f$ and @f$
330 * B_{Debye} @f$ are needed. The DebyeHuckel object uses two methods for
331 * specifying these quantities. The default method is to assume that @f$
332 * A_{Debye} @f$ is a constant, given in the initialization process, and stored
333 * in the member double, m_A_Debye. Optionally, a full water treatment may be
334 * employed that makes @f$ A_{Debye} @f$ a full function of *T* and *P*.
335 *
336 * @f[
337 * A_{Debye} = \frac{F e B_{Debye}}{8 \pi \epsilon R T} {\left( C_o \tilde{M}_o \right)}^{1/2}
338 * @f]
339 * where
340 * @f[
341 * B_{Debye} = \frac{F} {{(\frac{\epsilon R T}{2})}^{1/2}}
342 * @f]
343 * Therefore:
344 * @f[
345 * A_{Debye} = \frac{1}{8 \pi}
346 * {\left(\frac{2 N_a \rho_o}{1000}\right)}^{1/2}
347 * {\left(\frac{N_a e^2}{\epsilon R T }\right)}^{3/2}
348 * @f]
349 * where
350 * - @f$ N_a @f$ is Avogadro's number
351 * - @f$ \rho_w @f$ is the density of water
352 * - @f$ e @f$ is the electronic charge
353 * - @f$ \epsilon = K \epsilon_o @f$ is the permittivity of water
354 * - @f$ K @f$ is the dielectric constant of water
355 * - @f$ \epsilon_o @f$ is the permittivity of free space
356 * - @f$ \rho_o @f$ is the density of the solvent in its standard state.
357 *
358 * Nominal value at 298 K and 1 atm = 1.172576 (kg/gmol)^(1/2) based on:
359 * - @f$ \epsilon / \epsilon_0 @f$ = 78.54 (water at 25C)
360 * - T = 298.15 K
361 * - B_Debye = 3.28640E9 (kg/gmol)^(1/2) / m
362 *
363 * Currently, @f$ B_{Debye} @f$ is a constant in the model, specified either by
364 * a default water value, or through the input file. This may have to be looked
365 * at, in the future.
366 *
367 * Example phase and species definitions are given in the
368 * <a href="../../sphinx/html/yaml/phases.html#debye-huckel">YAML API Reference</a>.
369 *
370 * ## Application within Kinetics Managers
371 *
372 * For the time being, we have set the standard concentration for all species in
373 * this phase equal to the default concentration of the solvent at 298 K and 1
374 * atm. This means that the kinetics operator essentially works on an activities
375 * basis, with units specified as if it were on a concentration basis.
376 *
377 * For example, a bulk-phase binary reaction between liquid species j and k,
378 * producing a new liquid species l would have the following equation for its
379 * rate of progress variable, @f$ R^1 @f$, which has units of kmol m-3 s-1.
380 *
381 * @f[
382 * R^1 = k^1 C_j^a C_k^a = k^1 (C_o a_j) (C_o a_k)
383 * @f]
384 * where
385 * @f[
386 * C_j^a = C_o a_j \quad and \quad C_k^a = C_o a_k
387 * @f]
388 *
389 * @f$ C_j^a @f$ is the activity concentration of species j, and
390 * @f$ C_k^a @f$ is the activity concentration of species k. @f$ C_o @f$
391 * is the concentration of water at 298 K and 1 atm. @f$ a_j @f$ is the activity
392 * of species j at the current temperature and pressure and concentration of the
393 * liquid phase. @f$ k^1 @f$ has units of m3 kmol-1 s-1.
394 *
395 * The reverse rate constant can then be obtained from the law of microscopic
396 * reversibility and the equilibrium expression for the system.
397 *
398 * @f[
399 * \frac{a_j a_k}{ a_l} = K^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} )
400 * @f]
401 *
402 * @f$ K^{o,1} @f$ is the dimensionless form of the equilibrium constant.
403 *
404 * @f[
405 * R^{-1} = k^{-1} C_l^a = k^{-1} (C_o a_l)
406 * @f]
407 * where
408 * @f[
409 * k^{-1} = k^1 K^{o,1} C_o
410 * @f]
411 *
412 * @f$ k^{-1} @f$ has units of s-1.
413 */
415{
416public:
417 ~DebyeHuckel() override;
418
419 //! Full constructor for creating the phase.
420 /*!
421 * @param inputFile File name containing the definition of the phase.
422 * If blank, an empty phase will be created.
423 * @param id id attribute containing the name of the phase.
424 */
425 explicit DebyeHuckel(const string& inputFile="", const string& id="");
426
427 //! @name Utilities
428 //! @{
429
430 string type() const override {
431 return "Debye-Huckel";
432 }
433
434 //! @}
435 //! @name Mechanical Equation of State Properties
436 //!
437 //! In this equation of state implementation, the density is a function only
438 //! of the mole fractions. Therefore, it can't be an independent variable.
439 //! Instead, the pressure is used as the independent variable. Functions
440 //! which try to set the thermodynamic state by calling setDensity() will
441 //! cause an exception to be thrown.
442 //! @{
443
444protected:
445 void calcDensity() override;
446
447public:
448 //! @}
449 //! @name Activities, Standard States, and Activity Concentrations
450 //!
451 //! The activity @f$ a_k @f$ of a species in solution is related to the
452 //! chemical potential by @f[ \mu_k = \mu_k^0(T) + \hat R T \ln a_k. @f] The
453 //! quantity @f$ \mu_k^0(T,P) @f$ is the chemical potential at unit activity,
454 //! which depends only on temperature and the pressure. Activity is assumed
455 //! to be molality-based here.
456 //! @{
457
458 void getActivityConcentrations(double* c) const override;
459
460 //! Return the standard concentration for the kth species
461 /*!
462 * The standard concentration @f$ C^0_k @f$ used to normalize the activity
463 * (that is, generalized) concentration in kinetics calculations.
464 *
465 * For the time being, we will use the concentration of pure solvent for the
466 * the standard concentration of all species. This has the effect of making
467 * reaction rates based on the molality of species proportional to the
468 * molality of the species.
469 *
470 * @param k Optional parameter indicating the species. The default is to
471 * assume this refers to species 0.
472 * @return the standard Concentration in units of m^3/kmol
473 */
474 double standardConcentration(size_t k=0) const override;
475
476 //! Get the array of non-dimensional activities at the current solution
477 //! temperature, pressure, and solution concentration.
478 /*!
479 * (note solvent activity coefficient is on molar scale).
480 *
481 * @param ac Output vector of activities. Length: m_kk.
482 */
483 void getActivities(double* ac) const override;
484
485 //! Get the array of non-dimensional molality-based activity coefficients at
486 //! the current solution temperature, pressure, and solution concentration.
487 /*!
488 * note solvent is on molar scale. The solvent molar based activity
489 * coefficient is returned.
490 *
491 * Note, most of the work is done in an internal private routine
492 *
493 * @param acMolality Vector of Molality-based activity coefficients
494 * Length: m_kk
495 */
496 void getMolalityActivityCoefficients(double* acMolality) const override;
497
498 //! @}
499 //! @name Partial Molar Properties of the Solution
500 //! @{
501
502 //! Get the species chemical potentials. Units: J/kmol.
503 /*!
504 *
505 * This function returns a vector of chemical potentials of the species in
506 * solution.
507 *
508 * @f[
509 * \mu_k = \mu^{\triangle}_k(T,P) + R T \ln(\gamma_k^{\triangle} m_k)
510 * @f]
511 *
512 * @param mu Output vector of species chemical
513 * potentials. Length: m_kk. Units: J/kmol
514 */
515 void getChemPotentials(double* mu) const override;
516
517 //! Returns an array of partial molar enthalpies for the species
518 //! in the mixture. Units (J/kmol)
519 /*!
520 * For this phase, the partial molar enthalpies are equal to the
521 * standard state enthalpies modified by the derivative of the
522 * molality-based activity coefficient wrt temperature
523 *
524 * @f[
525 * \bar h_k(T,P) = h^{\triangle}_k(T,P) - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
526 * @f]
527 * The solvent partial molar enthalpy is equal to
528 * @f[
529 * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o}{dT}
530 * @f]
531 *
532 * The temperature dependence of the activity coefficients currently
533 * only occurs through the temperature dependence of the Debye constant.
534 *
535 * @param hbar Output vector of species partial molar enthalpies.
536 * Length: m_kk. units are J/kmol.
537 */
538 void getPartialMolarEnthalpies(double* hbar) const override;
539
540 //! Returns an array of partial molar entropies of the species in the
541 //! solution. Units: J/kmol/K.
542 /**
543 * Maxwell's equations provide an insight in how to calculate this
544 * (p.215 Smith and Van Ness)
545 * @f[
546 * \frac{d\mu_i}{dT} = -\bar{s}_i
547 * @f]
548 *
549 * For this phase, the partial molar entropies are equal to the SS species
550 * entropies plus the ideal solution contribution:
551 * @f[
552 * \bar s_k(T,P) = \hat s^0_k(T) - R \ln(M0 * molality[k])
553 * @f]
554 * @f[
555 * \bar s_{solvent}(T,P) = \hat s^0_{solvent}(T)
556 * - R ((xmolSolvent - 1.0) / xmolSolvent)
557 * @f]
558 *
559 * The reference-state pure-species entropies,@f$ \hat s^0_k(T) @f$, at the
560 * reference pressure, @f$ P_{ref} @f$, are computed by the species
561 * thermodynamic property manager. They are polynomial functions of
562 * temperature.
563 * @see MultiSpeciesThermo
564 *
565 * @param sbar Output vector of species partial molar entropies.
566 * Length = m_kk. units are J/kmol/K.
567 */
568 void getPartialMolarEntropies(double* sbar) const override;
569
570 void getPartialMolarCp(double* cpbar) const override;
571
572 //! Return an array of partial molar volumes for the species in the mixture.
573 //! Units: m^3/kmol.
574 /*!
575 * For this solution, the partial molar volumes are normally equal to the
576 * constant species molar volumes, except when the activity coefficients
577 * depend on pressure.
578 *
579 * The general relation is
580 *
581 * vbar_i = d(chemPot_i)/dP at const T, n
582 * = V0_i + d(Gex)/dP)_T,M
583 * = V0_i + RT d(lnActCoeffi)dP _T,M
584 *
585 * @param vbar Output vector of species partial molar volumes.
586 * Length = m_kk. units are m^3/kmol.
587 */
588 void getPartialMolarVolumes(double* vbar) const override;
589
590 //! @}
591
592 /*
593 * -------------- Utilities -------------------------------
594 */
595
596 bool addSpecies(shared_ptr<Species> spec) override;
597 void initThermo() override;
598 void getParameters(AnyMap& phaseNode) const override;
599 void getSpeciesParameters(const string& name, AnyMap& speciesNode) const override;
600
601 //! Return the Debye Huckel constant as a function of temperature
602 //! and pressure (Units = sqrt(kg/gmol))
603 /*!
604 * The default is to assume that it is constant, given in the
605 * initialization process, and stored in the member double, m_A_Debye.
606 * Optionally, a full water treatment may be employed that makes
607 * @f$ A_{Debye} @f$ a full function of T and P.
608 *
609 * @f[
610 * A_{Debye} = \frac{F e B_{Debye}}{8 \pi \epsilon R T} {\left( C_o \tilde{M}_o \right)}^{1/2}
611 * @f]
612 * where
613 * @f[
614 * B_{Debye} = \frac{F} {{(\frac{\epsilon R T}{2})}^{1/2}}
615 * @f]
616 * Therefore:
617 * @f[
618 * A_{Debye} = \frac{1}{8 \pi}
619 * {\left(\frac{2 N_a \rho_o}{1000}\right)}^{1/2}
620 * {\left(\frac{N_a e^2}{\epsilon R T }\right)}^{3/2}
621 * @f]
622 *
623 * where
624 * - Units = sqrt(kg/gmol)
625 * - @f$ N_a @f$ is Avogadro's number
626 * - @f$ \rho_w @f$ is the density of water
627 * - @f$ e @f$ is the electronic charge
628 * - @f$ \epsilon = K \epsilon_o @f$ is the permittivity of water
629 * - @f$ K @f$ is the dielectric constant of water,
630 * - @f$ \epsilon_o @f$ is the permittivity of free space.
631 * - @f$ \rho_o @f$ is the density of the solvent in its standard state.
632 *
633 * Nominal value at 298 K and 1 atm = 1.172576 (kg/gmol)^(1/2)
634 * based on:
635 * - @f$ \epsilon / \epsilon_0 @f$ = 78.54 (water at 25C)
636 * - T = 298.15 K
637 * - B_Debye = 3.28640E9 (kg/gmol)^(1/2)/m
638 *
639 * @param temperature Temperature in kelvin. Defaults to -1, in which
640 * case the temperature of the phase is assumed.
641 * @param pressure Pressure (Pa). Defaults to -1, in which
642 * case the pressure of the phase is assumed.
643 */
644 virtual double A_Debye_TP(double temperature = -1.0,
645 double pressure = -1.0) const;
646
647 //! Value of the derivative of the Debye Huckel constant with
648 //! respect to temperature.
649 /*!
650 * This is a function of temperature and pressure. See A_Debye_TP() for
651 * a definition of @f$ A_{Debye} @f$.
652 *
653 * Units = sqrt(kg/gmol) K-1
654 *
655 * @param temperature Temperature in kelvin. Defaults to -1, in which
656 * case the temperature of the phase is assumed.
657 * @param pressure Pressure (Pa). Defaults to -1, in which
658 * case the pressure of the phase is assumed.
659 */
660 virtual double dA_DebyedT_TP(double temperature = -1.0,
661 double pressure = -1.0) const;
662
663 //! Value of the 2nd derivative of the Debye Huckel constant with
664 //! respect to temperature as a function of temperature and pressure.
665 /*!
666 * This is a function of temperature and pressure. See A_Debye_TP() for
667 * a definition of @f$ A_{Debye} @f$.
668 *
669 * Units = sqrt(kg/gmol) K-2
670 *
671 * @param temperature Temperature in kelvin. Defaults to -1, in which
672 * case the temperature of the phase is assumed.
673 * @param pressure Pressure (Pa). Defaults to -1, in which
674 * case the pressure of the phase is assumed.
675 */
676 virtual double d2A_DebyedT2_TP(double temperature = -1.0,
677 double pressure = -1.0) const;
678
679 //! Value of the derivative of the Debye Huckel constant with
680 //! respect to pressure, as a function of temperature and pressure.
681 /*!
682 * This is a function of temperature and pressure. See A_Debye_TP() for
683 * a definition of @f$ A_{Debye} @f$.
684 *
685 * Units = sqrt(kg/gmol) Pa-1
686 *
687 * @param temperature Temperature in kelvin. Defaults to -1, in which
688 * case the temperature of the phase is assumed.
689 * @param pressure Pressure (Pa). Defaults to -1, in which
690 * case the pressure of the phase is assumed.
691 */
692 virtual double dA_DebyedP_TP(double temperature = -1.0,
693 double pressure = -1.0) const;
694
695 //! Reports the ionic radius of the kth species
696 /*!
697 * @param k species index.
698 */
699 double AionicRadius(int k = 0) const;
700
701 //! Set the DebyeHuckel parameterization form. Must be one of
702 //! 'dilute-limit', 'B-dot-with-variable-a', 'B-dot-with-common-a',
703 //! 'beta_ij', or 'Pitzer-with-beta_ij'.
704 void setDebyeHuckelModel(const string& form);
705
706 //! Returns the form of the Debye-Huckel parameterization used
707 int formDH() const {
708 return m_formDH;
709 }
710
711 //! Set the A_Debye parameter. If a negative value is provided, enables
712 //! calculation of A_Debye using the detailed water equation of state.
713 void setA_Debye(double A);
714
715 void setB_Debye(double B) { m_B_Debye = B; }
716 void setB_dot(double bdot);
717 void setMaxIonicStrength(double Imax) { m_maxIionicStrength = Imax; }
718 void useHelgesonFixedForm(bool mode=true) { m_useHelgesonFixedForm = mode; }
719
720 //! Set the default ionic radius [m] for each species
721 void setDefaultIonicRadius(double value);
722
723 //! Set the value for the beta interaction between species sp1 and sp2.
724 void setBeta(const string& sp1, const string& sp2, double value);
725
726 //! Returns a reference to M_Beta_ij
728 return m_Beta_ij;
729 }
730
731private:
732 //! Static function that implements the non-polar species salt-out
733 //! modifications.
734 /*!
735 * Returns the calculated activity coefficients.
736 *
737 * @param IionicMolality Value of the ionic molality (sqrt(gmol/kg))
738 */
739 static double _nonpolarActCoeff(double IionicMolality);
740
741 //! Formula for the osmotic coefficient that occurs in the GWB.
742 /*!
743 * It is originally from Helgeson for a variable NaCl brine. It's to be
744 * used with extreme caution.
745 */
746 double _osmoticCoeffHelgesonFixedForm() const;
747
748 //! Formula for the log of the water activity that occurs in the GWB.
749 /*!
750 * It is originally from Helgeson for a variable NaCl brine. It's to be
751 * used with extreme caution.
752 */
754
755protected:
756 //! form of the Debye-Huckel parameterization used in the model.
757 /*!
758 * The options are described at the top of this document,
759 * and in the general documentation.
760 * The list is repeated here:
761 *
762 * DHFORM_DILUTE_LIMIT = 0 (default)
763 * DHFORM_BDOT_AK = 1
764 * DHFORM_BDOT_AUNIFORM = 2
765 * DHFORM_BETAIJ = 3
766 * DHFORM_PITZER_BETAIJ = 4
767 */
768 int m_formDH = DHFORM_DILUTE_LIMIT;
769
770 //! Vector containing the electrolyte species type
771 /*!
772 * The possible types are:
773 * - solvent
774 * - Charged Species
775 * - weakAcidAssociated
776 * - strongAcidAssociated
777 * - polarNeutral
778 * - nonpolarNeutral
779 * .
780 */
782
783 //! a_k = Size of the ionic species in the DH formulation. units = meters
784 vector<double> m_Aionic;
785
786 //! Default ionic radius for species where it is not specified
787 double m_Aionic_default = NAN;
788
789 //! Current value of the ionic strength on the molality scale
790 mutable double m_IionicMolality = 0.0;
791
792 //! Maximum value of the ionic strength allowed in the calculation of the
793 //! activity coefficients.
795
796public:
797 //! If true, then the fixed for of Helgeson's activity for water is used
798 //! instead of the rigorous form obtained from Gibbs-Duhem relation. This
799 //! should be used with caution, and is really only included as a validation
800 //! exercise.
802protected:
803 //! Stoichiometric ionic strength on the molality scale
804 mutable double m_IionicMolalityStoich = 0.0;
805
806public:
807 /**
808 * Form of the constant outside the Debye-Huckel term
809 * called A. It's normally a function of temperature
810 * and pressure. However, it can be set from the
811 * input file in order to aid in numerical comparisons.
812 * Acceptable forms:
813 *
814 * A_DEBYE_CONST 0
815 * A_DEBYE_WATER 1
816 *
817 * The A_DEBYE_WATER form may be used for water solvents
818 * with needs to cover varying temperatures and pressures.
819 * Note, the dielectric constant of water is a relatively
820 * strong function of T, and its variability must be
821 * accounted for,
822 */
823 int m_form_A_Debye = A_DEBYE_CONST;
824
825protected:
826 //! Current value of the Debye Constant, A_Debye
827 /**
828 * A_Debye -> this expression appears on the top of the ln actCoeff term in
829 * the general Debye-Huckel expression It depends on temperature
830 * and pressure.
831 *
832 * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
833 *
834 * Units = sqrt(kg/gmol)
835 *
836 * Nominal value(298K, atm) = 1.172576 sqrt(kg/gmol)
837 * based on:
838 * epsilon/epsilon_0 = 78.54
839 * (water at 25C)
840 * T = 298.15 K
841 * B_Debye = 3.28640E9 sqrt(kg/gmol)/m
842 *
843 * note in Pitzer's nomenclature, A_phi = A_Debye/3.0
844 */
845 mutable double m_A_Debye;
846
847 //! Current value of the constant that appears in the denominator
848 /**
849 * B_Debye -> this expression appears on the bottom of the ln actCoeff term
850 * in the general Debye-Huckel expression It depends on
851 * temperature
852 *
853 * B_Bebye = F / sqrt( epsilon R T / 2 )
854 *
855 * Units = sqrt(kg/gmol) / m
856 *
857 * Nominal value = 3.28640E9 sqrt(kg/gmol) / m
858 * based on:
859 * epsilon/epsilon_0 = 78.54
860 * (water at 25C)
861 * T = 298.15 K
862 */
863 double m_B_Debye;
864
865 //! Array of B_Dot values
866 /**
867 * This expression is an extension of the Debye-Huckel expression used
868 * in some formulations to extend DH to higher molalities. B_dot is
869 * specific to the major ionic pair.
870 */
871 vector<double> m_B_Dot;
872
873 //! Pointer to the Water standard state object
874 /*!
875 * derived from the equation of state for water.
876 */
878
879 //! Storage for the density of water's standard state
880 /*!
881 * Density depends on temperature and pressure.
882 */
883 double m_densWaterSS = 1000.0;
884
885 //! Pointer to the water property calculator
886 unique_ptr<WaterProps> m_waterProps;
887
888 /**
889 * Stoichiometric species charge -> This is for calculations
890 * of the ionic strength which ignore ion-ion pairing into
891 * neutral molecules. The Stoichiometric species charge is the
892 * charge of one of the ion that would occur if the species broke
893 * into two charged ion pairs.
894 * NaCl -> m_speciesCharge_Stoich = -1;
895 * HSO4- -> H+ + SO42- = -2
896 * -> The other charge is calculated.
897 * For species that aren't ion pairs, it's equal to the
898 * m_speciesCharge[] value.
899 */
901
902 /**
903 * Array of 2D data used in the DHFORM_BETAIJ formulation
904 * Beta_ij.value(i,j) is the coefficient of the jth species
905 * for the specification of the chemical potential of the ith
906 * species.
907 */
909
910 //! Logarithm of the activity coefficients on the molality scale.
911 /*!
912 * mutable because we change this if the composition or temperature or
913 * pressure changes.
914 */
915 mutable vector<double> m_lnActCoeffMolal;
916
917 //! Derivative of log act coeff wrt T
918 mutable vector<double> m_dlnActCoeffMolaldT;
919
920 //! 2nd Derivative of log act coeff wrt T
921 mutable vector<double> m_d2lnActCoeffMolaldT2;
922
923 //! Derivative of log act coeff wrt P
924 mutable vector<double> m_dlnActCoeffMolaldP;
925
926private:
927 //! Calculate the log activity coefficients
928 /*!
929 * This function updates the internally stored natural logarithm of the
930 * molality activity coefficients. This is the main routine for
931 * implementing the activity coefficient formulation.
932 */
933 void s_update_lnMolalityActCoeff() const;
934
935 //! Calculation of temperature derivative of activity coefficient
936 /*!
937 * Using internally stored values, this function calculates the temperature
938 * derivative of the logarithm of the activity coefficient for all species
939 * in the mechanism.
940 *
941 * We assume that the activity coefficients are current in this routine. The
942 * solvent activity coefficient is on the molality scale. Its derivative is
943 * too.
944 */
946
947 //! Calculate the temperature 2nd derivative of the activity coefficient
948 /*!
949 * Using internally stored values, this function calculates the temperature
950 * 2nd derivative of the logarithm of the activity coefficient for all
951 * species in the mechanism.
952 *
953 * We assume that the activity coefficients are current in this routine.
954 * Solvent activity coefficient is on the molality scale. Its derivatives
955 * are too.
956 */
958
959 //! Calculate the pressure derivative of the activity coefficient
960 /*!
961 * Using internally stored values, this function calculates the pressure
962 * derivative of the logarithm of the activity coefficient for all species
963 * in the mechanism.
964 *
965 * We assume that the activity coefficients, molalities, and A_Debye are
966 * current. Solvent activity coefficient is on the molality scale. Its
967 * derivatives are too.
968 */
970};
971
972}
973
974#endif
Header file for class Cantera::Array2D.
Header for intermediate ThermoPhase object for phases which employ molality based activity coefficien...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:432
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
Class DebyeHuckel represents a dilute liquid electrolyte phase which obeys the Debye Huckel formulati...
Array2D m_Beta_ij
Array of 2D data used in the DHFORM_BETAIJ formulation Beta_ij.value(i,j) is the coefficient of the j...
Array2D & get_Beta_ij()
Returns a reference to M_Beta_ij.
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
int m_formDH
form of the Debye-Huckel parameterization used in the model.
double m_A_Debye
Current value of the Debye Constant, A_Debye.
virtual double d2A_DebyedT2_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the 2nd derivative of the Debye Huckel constant with respect to temperature as a function of...
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
vector< double > m_B_Dot
Array of B_Dot values.
double m_IionicMolality
Current value of the ionic strength on the molality scale.
double _osmoticCoeffHelgesonFixedForm() const
Formula for the osmotic coefficient that occurs in the GWB.
double m_densWaterSS
Storage for the density of water's standard state.
void s_update_d2lnMolalityActCoeff_dT2() const
Calculate the temperature 2nd derivative of the activity coefficient.
void getSpeciesParameters(const string &name, AnyMap &speciesNode) const override
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
bool m_useHelgesonFixedForm
If true, then the fixed for of Helgeson's activity for water is used instead of the rigorous form obt...
static double _nonpolarActCoeff(double IionicMolality)
Static function that implements the non-polar species salt-out modifications.
vector< int > m_electrolyteSpeciesType
Vector containing the electrolyte species type.
int formDH() const
Returns the form of the Debye-Huckel parameterization used.
double m_B_Debye
Current value of the constant that appears in the denominator.
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
string type() const override
String indicating the thermodynamic model implemented.
double AionicRadius(int k=0) const
Reports the ionic radius of the kth species.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void getActivityConcentrations(double *c) const override
This method returns an array of generalized concentrations.
void s_update_dlnMolalityActCoeff_dT() const
Calculation of temperature derivative of activity coefficient.
void s_update_lnMolalityActCoeff() const
Calculate the log activity coefficients.
vector< double > m_Aionic
a_k = Size of the ionic species in the DH formulation. units = meters
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
void setA_Debye(double A)
Set the A_Debye parameter.
void calcDensity() override
Calculate the density of the mixture using the partial molar volumes and mole fractions as input.
vector< double > m_dlnActCoeffMolaldT
Derivative of log act coeff wrt T.
double m_Aionic_default
Default ionic radius for species where it is not specified.
void setDebyeHuckelModel(const string &form)
Set the DebyeHuckel parameterization form.
vector< double > m_speciesCharge_Stoich
Stoichiometric species charge -> This is for calculations of the ionic strength which ignore ion-ion ...
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
double _lnactivityWaterHelgesonFixedForm() const
Formula for the log of the water activity that occurs in the GWB.
void setBeta(const string &sp1, const string &sp2, double value)
Set the value for the beta interaction between species sp1 and sp2.
void getActivities(double *ac) const override
Get the array of non-dimensional activities at the current solution temperature, pressure,...
virtual double dA_DebyedT_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to temperature.
void getPartialMolarCp(double *cpbar) const override
Return an array of partial molar heat capacities for the species in the mixture.
void setDefaultIonicRadius(double value)
Set the default ionic radius [m] for each species.
vector< double > m_d2lnActCoeffMolaldT2
2nd Derivative of log act coeff wrt T
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
void s_update_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the activity coefficient.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
PDSS_Water * m_waterSS
Pointer to the Water standard state object.
double m_IionicMolalityStoich
Stoichiometric ionic strength on the molality scale.
vector< double > m_lnActCoeffMolal
Logarithm of the activity coefficients on the molality scale.
vector< double > m_dlnActCoeffMolaldP
Derivative of log act coeff wrt P.
void getMolalityActivityCoefficients(double *acMolality) const override
Get the array of non-dimensional molality-based activity coefficients at the current solution tempera...
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
virtual double A_Debye_TP(double temperature=-1.0, double pressure=-1.0) const
Return the Debye Huckel constant as a function of temperature and pressure (Units = sqrt(kg/gmol))
virtual double dA_DebyedP_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to pressure, as a function of tempe...
MolalityVPSSTP is a derived class of ThermoPhase that handles variable pressure standard state method...
Class for the liquid water pressure dependent standard state.
Definition PDSS_Water.h:50
double temperature() const
Temperature (K).
Definition Phase.h:563
string name() const
Return the name of the phase.
Definition Phase.cpp:20
double pressure() const override
Returns the current pressure of the phase.
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595