Cantera  3.1.0a1
Loading...
Searching...
No Matches
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 Molar Thermodynamic Properties of the Solution
436 //! @{
437
438 double enthalpy_mole() const override;
439
440 //! Molar entropy. Units: J/kmol/K.
441 /**
442 * For an ideal, constant partial molar volume solution mixture with
443 * pure species phases which exhibit zero volume expansivity:
444 * @f[
445 * \hat s(T, P, X_k) = \sum_k X_k \hat s^0_k(T)
446 * - \hat R \sum_k X_k \ln(X_k)
447 * @f]
448 * The reference-state pure-species entropies
449 * @f$ \hat s^0_k(T,p_{ref}) @f$ are computed by the
450 * species thermodynamic
451 * property manager. The pure species entropies are independent of
452 * temperature since the volume expansivities are equal to zero.
453 * @see MultiSpeciesThermo
454 */
455 double entropy_mole() const override;
456
457 double gibbs_mole() const override;
458 double cp_mole() const override;
459
460 //! @}
461 //! @name Mechanical Equation of State Properties
462 //!
463 //! In this equation of state implementation, the density is a function only
464 //! of the mole fractions. Therefore, it can't be an independent variable.
465 //! Instead, the pressure is used as the independent variable. Functions
466 //! which try to set the thermodynamic state by calling setDensity() will
467 //! cause an exception to be thrown.
468 //! @{
469
470protected:
471 void calcDensity() override;
472
473public:
474 //! @}
475 //! @name Activities, Standard States, and Activity Concentrations
476 //!
477 //! The activity @f$ a_k @f$ of a species in solution is related to the
478 //! chemical potential by @f[ \mu_k = \mu_k^0(T) + \hat R T \ln a_k. @f] The
479 //! quantity @f$ \mu_k^0(T,P) @f$ is the chemical potential at unit activity,
480 //! which depends only on temperature and the pressure. Activity is assumed
481 //! to be molality-based here.
482 //! @{
483
484 void getActivityConcentrations(double* c) const override;
485
486 //! Return the standard concentration for the kth species
487 /*!
488 * The standard concentration @f$ C^0_k @f$ used to normalize the activity
489 * (that is, generalized) concentration in kinetics calculations.
490 *
491 * For the time being, we will use the concentration of pure solvent for the
492 * the standard concentration of all species. This has the effect of making
493 * reaction rates based on the molality of species proportional to the
494 * molality of the species.
495 *
496 * @param k Optional parameter indicating the species. The default is to
497 * assume this refers to species 0.
498 * @return the standard Concentration in units of m^3/kmol
499 */
500 double standardConcentration(size_t k=0) const override;
501
502 //! Get the array of non-dimensional activities at the current solution
503 //! temperature, pressure, and solution concentration.
504 /*!
505 * (note solvent activity coefficient is on molar scale).
506 *
507 * @param ac Output vector of activities. Length: m_kk.
508 */
509 void getActivities(double* ac) const override;
510
511 //! Get the array of non-dimensional molality-based activity coefficients at
512 //! the current solution temperature, pressure, and solution concentration.
513 /*!
514 * note solvent is on molar scale. The solvent molar based activity
515 * coefficient is returned.
516 *
517 * Note, most of the work is done in an internal private routine
518 *
519 * @param acMolality Vector of Molality-based activity coefficients
520 * Length: m_kk
521 */
522 void getMolalityActivityCoefficients(double* acMolality) const override;
523
524 //! @}
525 //! @name Partial Molar Properties of the Solution
526 //! @{
527
528 //! Get the species chemical potentials. Units: J/kmol.
529 /*!
530 *
531 * This function returns a vector of chemical potentials of the species in
532 * solution.
533 *
534 * @f[
535 * \mu_k = \mu^{\triangle}_k(T,P) + R T \ln(\gamma_k^{\triangle} m_k)
536 * @f]
537 *
538 * @param mu Output vector of species chemical
539 * potentials. Length: m_kk. Units: J/kmol
540 */
541 void getChemPotentials(double* mu) const override;
542
543 //! Returns an array of partial molar enthalpies for the species
544 //! in the mixture. Units (J/kmol)
545 /*!
546 * For this phase, the partial molar enthalpies are equal to the
547 * standard state enthalpies modified by the derivative of the
548 * molality-based activity coefficient wrt temperature
549 *
550 * @f[
551 * \bar h_k(T,P) = h^{\triangle}_k(T,P) - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
552 * @f]
553 * The solvent partial molar enthalpy is equal to
554 * @f[
555 * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o}{dT}
556 * @f]
557 *
558 * The temperature dependence of the activity coefficients currently
559 * only occurs through the temperature dependence of the Debye constant.
560 *
561 * @param hbar Output vector of species partial molar enthalpies.
562 * Length: m_kk. units are J/kmol.
563 */
564 void getPartialMolarEnthalpies(double* hbar) const override;
565
566 //! Returns an array of partial molar entropies of the species in the
567 //! solution. Units: J/kmol/K.
568 /**
569 * Maxwell's equations provide an insight in how to calculate this
570 * (p.215 Smith and Van Ness)
571 * @f[
572 * \frac{d\mu_i}{dT} = -\bar{s}_i
573 * @f]
574 *
575 * For this phase, the partial molar entropies are equal to the SS species
576 * entropies plus the ideal solution contribution:
577 * @f[
578 * \bar s_k(T,P) = \hat s^0_k(T) - R \ln(M0 * molality[k])
579 * @f]
580 * @f[
581 * \bar s_{solvent}(T,P) = \hat s^0_{solvent}(T)
582 * - R ((xmolSolvent - 1.0) / xmolSolvent)
583 * @f]
584 *
585 * The reference-state pure-species entropies,@f$ \hat s^0_k(T) @f$, at the
586 * reference pressure, @f$ P_{ref} @f$, are computed by the species
587 * thermodynamic property manager. They are polynomial functions of
588 * temperature.
589 * @see MultiSpeciesThermo
590 *
591 * @param sbar Output vector of species partial molar entropies.
592 * Length = m_kk. units are J/kmol/K.
593 */
594 void getPartialMolarEntropies(double* sbar) const override;
595
596 void getPartialMolarCp(double* cpbar) const override;
597
598 //! Return an array of partial molar volumes for the species in the mixture.
599 //! Units: m^3/kmol.
600 /*!
601 * For this solution, the partial molar volumes are normally equal to the
602 * constant species molar volumes, except when the activity coefficients
603 * depend on pressure.
604 *
605 * The general relation is
606 *
607 * vbar_i = d(chemPot_i)/dP at const T, n
608 * = V0_i + d(Gex)/dP)_T,M
609 * = V0_i + RT d(lnActCoeffi)dP _T,M
610 *
611 * @param vbar Output vector of species partial molar volumes.
612 * Length = m_kk. units are m^3/kmol.
613 */
614 void getPartialMolarVolumes(double* vbar) const override;
615
616 //! @}
617
618 /*
619 * -------------- Utilities -------------------------------
620 */
621
622 bool addSpecies(shared_ptr<Species> spec) override;
623 void initThermo() override;
624 void getParameters(AnyMap& phaseNode) const override;
625 void getSpeciesParameters(const string& name, AnyMap& speciesNode) const override;
626
627 //! Return the Debye Huckel constant as a function of temperature
628 //! and pressure (Units = sqrt(kg/gmol))
629 /*!
630 * The default is to assume that it is constant, given in the
631 * initialization process, and stored in the member double, m_A_Debye.
632 * Optionally, a full water treatment may be employed that makes
633 * @f$ A_{Debye} @f$ a full function of T and P.
634 *
635 * @f[
636 * A_{Debye} = \frac{F e B_{Debye}}{8 \pi \epsilon R T} {\left( C_o \tilde{M}_o \right)}^{1/2}
637 * @f]
638 * where
639 * @f[
640 * B_{Debye} = \frac{F} {{(\frac{\epsilon R T}{2})}^{1/2}}
641 * @f]
642 * Therefore:
643 * @f[
644 * A_{Debye} = \frac{1}{8 \pi}
645 * {\left(\frac{2 N_a \rho_o}{1000}\right)}^{1/2}
646 * {\left(\frac{N_a e^2}{\epsilon R T }\right)}^{3/2}
647 * @f]
648 *
649 * where
650 * - Units = sqrt(kg/gmol)
651 * - @f$ N_a @f$ is Avogadro's number
652 * - @f$ \rho_w @f$ is the density of water
653 * - @f$ e @f$ is the electronic charge
654 * - @f$ \epsilon = K \epsilon_o @f$ is the permittivity of water
655 * - @f$ K @f$ is the dielectric constant of water,
656 * - @f$ \epsilon_o @f$ is the permittivity of free space.
657 * - @f$ \rho_o @f$ is the density of the solvent in its standard state.
658 *
659 * Nominal value at 298 K and 1 atm = 1.172576 (kg/gmol)^(1/2)
660 * based on:
661 * - @f$ \epsilon / \epsilon_0 @f$ = 78.54 (water at 25C)
662 * - T = 298.15 K
663 * - B_Debye = 3.28640E9 (kg/gmol)^(1/2)/m
664 *
665 * @param temperature Temperature in kelvin. Defaults to -1, in which
666 * case the temperature of the phase is assumed.
667 * @param pressure Pressure (Pa). Defaults to -1, in which
668 * case the pressure of the phase is assumed.
669 */
670 virtual double A_Debye_TP(double temperature = -1.0,
671 double pressure = -1.0) const;
672
673 //! Value of the derivative of the Debye Huckel constant with
674 //! respect to temperature.
675 /*!
676 * This is a function of temperature and pressure. See A_Debye_TP() for
677 * a definition of @f$ A_{Debye} @f$.
678 *
679 * Units = sqrt(kg/gmol) K-1
680 *
681 * @param temperature Temperature in kelvin. Defaults to -1, in which
682 * case the temperature of the phase is assumed.
683 * @param pressure Pressure (Pa). Defaults to -1, in which
684 * case the pressure of the phase is assumed.
685 */
686 virtual double dA_DebyedT_TP(double temperature = -1.0,
687 double pressure = -1.0) const;
688
689 //! Value of the 2nd derivative of the Debye Huckel constant with
690 //! respect to temperature as a function of temperature and pressure.
691 /*!
692 * This is a function of temperature and pressure. See A_Debye_TP() for
693 * a definition of @f$ A_{Debye} @f$.
694 *
695 * Units = sqrt(kg/gmol) K-2
696 *
697 * @param temperature Temperature in kelvin. Defaults to -1, in which
698 * case the temperature of the phase is assumed.
699 * @param pressure Pressure (Pa). Defaults to -1, in which
700 * case the pressure of the phase is assumed.
701 */
702 virtual double d2A_DebyedT2_TP(double temperature = -1.0,
703 double pressure = -1.0) const;
704
705 //! Value of the derivative of the Debye Huckel constant with
706 //! respect to pressure, as a function of temperature and pressure.
707 /*!
708 * This is a function of temperature and pressure. See A_Debye_TP() for
709 * a definition of @f$ A_{Debye} @f$.
710 *
711 * Units = sqrt(kg/gmol) Pa-1
712 *
713 * @param temperature Temperature in kelvin. Defaults to -1, in which
714 * case the temperature of the phase is assumed.
715 * @param pressure Pressure (Pa). Defaults to -1, in which
716 * case the pressure of the phase is assumed.
717 */
718 virtual double dA_DebyedP_TP(double temperature = -1.0,
719 double pressure = -1.0) const;
720
721 //! Reports the ionic radius of the kth species
722 /*!
723 * @param k species index.
724 */
725 double AionicRadius(int k = 0) const;
726
727 //! Set the DebyeHuckel parameterization form. Must be one of
728 //! 'dilute-limit', 'B-dot-with-variable-a', 'B-dot-with-common-a',
729 //! 'beta_ij', or 'Pitzer-with-beta_ij'.
730 void setDebyeHuckelModel(const string& form);
731
732 //! Returns the form of the Debye-Huckel parameterization used
733 int formDH() const {
734 return m_formDH;
735 }
736
737 //! Set the A_Debye parameter. If a negative value is provided, enables
738 //! calculation of A_Debye using the detailed water equation of state.
739 void setA_Debye(double A);
740
741 void setB_Debye(double B) { m_B_Debye = B; }
742 void setB_dot(double bdot);
743 void setMaxIonicStrength(double Imax) { m_maxIionicStrength = Imax; }
744 void useHelgesonFixedForm(bool mode=true) { m_useHelgesonFixedForm = mode; }
745
746 //! Set the default ionic radius [m] for each species
747 void setDefaultIonicRadius(double value);
748
749 //! Set the value for the beta interaction between species sp1 and sp2.
750 void setBeta(const string& sp1, const string& sp2, double value);
751
752 //! Returns a reference to M_Beta_ij
754 return m_Beta_ij;
755 }
756
757private:
758 //! Static function that implements the non-polar species salt-out
759 //! modifications.
760 /*!
761 * Returns the calculated activity coefficients.
762 *
763 * @param IionicMolality Value of the ionic molality (sqrt(gmol/kg))
764 */
765 static double _nonpolarActCoeff(double IionicMolality);
766
767 //! Formula for the osmotic coefficient that occurs in the GWB.
768 /*!
769 * It is originally from Helgeson for a variable NaCl brine. It's to be
770 * used with extreme caution.
771 */
772 double _osmoticCoeffHelgesonFixedForm() const;
773
774 //! Formula for the log of the water activity that occurs in the GWB.
775 /*!
776 * It is originally from Helgeson for a variable NaCl brine. It's to be
777 * used with extreme caution.
778 */
780
781protected:
782 //! form of the Debye-Huckel parameterization used in the model.
783 /*!
784 * The options are described at the top of this document,
785 * and in the general documentation.
786 * The list is repeated here:
787 *
788 * DHFORM_DILUTE_LIMIT = 0 (default)
789 * DHFORM_BDOT_AK = 1
790 * DHFORM_BDOT_AUNIFORM = 2
791 * DHFORM_BETAIJ = 3
792 * DHFORM_PITZER_BETAIJ = 4
793 */
794 int m_formDH = DHFORM_DILUTE_LIMIT;
795
796 //! Vector containing the electrolyte species type
797 /*!
798 * The possible types are:
799 * - solvent
800 * - Charged Species
801 * - weakAcidAssociated
802 * - strongAcidAssociated
803 * - polarNeutral
804 * - nonpolarNeutral
805 * .
806 */
808
809 //! a_k = Size of the ionic species in the DH formulation. units = meters
810 vector<double> m_Aionic;
811
812 //! Default ionic radius for species where it is not specified
813 double m_Aionic_default = NAN;
814
815 //! Current value of the ionic strength on the molality scale
816 mutable double m_IionicMolality = 0.0;
817
818 //! Maximum value of the ionic strength allowed in the calculation of the
819 //! activity coefficients.
821
822public:
823 //! If true, then the fixed for of Helgeson's activity for water is used
824 //! instead of the rigorous form obtained from Gibbs-Duhem relation. This
825 //! should be used with caution, and is really only included as a validation
826 //! exercise.
828protected:
829 //! Stoichiometric ionic strength on the molality scale
830 mutable double m_IionicMolalityStoich = 0.0;
831
832public:
833 /**
834 * Form of the constant outside the Debye-Huckel term
835 * called A. It's normally a function of temperature
836 * and pressure. However, it can be set from the
837 * input file in order to aid in numerical comparisons.
838 * Acceptable forms:
839 *
840 * A_DEBYE_CONST 0
841 * A_DEBYE_WATER 1
842 *
843 * The A_DEBYE_WATER form may be used for water solvents
844 * with needs to cover varying temperatures and pressures.
845 * Note, the dielectric constant of water is a relatively
846 * strong function of T, and its variability must be
847 * accounted for,
848 */
849 int m_form_A_Debye = A_DEBYE_CONST;
850
851protected:
852 //! Current value of the Debye Constant, A_Debye
853 /**
854 * A_Debye -> this expression appears on the top of the ln actCoeff term in
855 * the general Debye-Huckel expression It depends on temperature
856 * and pressure.
857 *
858 * A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
859 *
860 * Units = sqrt(kg/gmol)
861 *
862 * Nominal value(298K, atm) = 1.172576 sqrt(kg/gmol)
863 * based on:
864 * epsilon/epsilon_0 = 78.54
865 * (water at 25C)
866 * T = 298.15 K
867 * B_Debye = 3.28640E9 sqrt(kg/gmol)/m
868 *
869 * note in Pitzer's nomenclature, A_phi = A_Debye/3.0
870 */
871 mutable double m_A_Debye;
872
873 //! Current value of the constant that appears in the denominator
874 /**
875 * B_Debye -> this expression appears on the bottom of the ln actCoeff term
876 * in the general Debye-Huckel expression It depends on
877 * temperature
878 *
879 * B_Bebye = F / sqrt( epsilon R T / 2 )
880 *
881 * Units = sqrt(kg/gmol) / m
882 *
883 * Nominal value = 3.28640E9 sqrt(kg/gmol) / m
884 * based on:
885 * epsilon/epsilon_0 = 78.54
886 * (water at 25C)
887 * T = 298.15 K
888 */
889 double m_B_Debye;
890
891 //! Array of B_Dot values
892 /**
893 * This expression is an extension of the Debye-Huckel expression used
894 * in some formulations to extend DH to higher molalities. B_dot is
895 * specific to the major ionic pair.
896 */
897 vector<double> m_B_Dot;
898
899 //! Pointer to the Water standard state object
900 /*!
901 * derived from the equation of state for water.
902 */
904
905 //! Storage for the density of water's standard state
906 /*!
907 * Density depends on temperature and pressure.
908 */
909 double m_densWaterSS = 1000.0;
910
911 //! Pointer to the water property calculator
912 unique_ptr<WaterProps> m_waterProps;
913
914 //! vector of size m_kk, used as a temporary holding area.
915 mutable vector<double> m_tmpV;
916
917 /**
918 * Stoichiometric species charge -> This is for calculations
919 * of the ionic strength which ignore ion-ion pairing into
920 * neutral molecules. The Stoichiometric species charge is the
921 * charge of one of the ion that would occur if the species broke
922 * into two charged ion pairs.
923 * NaCl -> m_speciesCharge_Stoich = -1;
924 * HSO4- -> H+ + SO42- = -2
925 * -> The other charge is calculated.
926 * For species that aren't ion pairs, it's equal to the
927 * m_speciesCharge[] value.
928 */
930
931 /**
932 * Array of 2D data used in the DHFORM_BETAIJ formulation
933 * Beta_ij.value(i,j) is the coefficient of the jth species
934 * for the specification of the chemical potential of the ith
935 * species.
936 */
938
939 //! Logarithm of the activity coefficients on the molality scale.
940 /*!
941 * mutable because we change this if the composition or temperature or
942 * pressure changes.
943 */
944 mutable vector<double> m_lnActCoeffMolal;
945
946 //! Derivative of log act coeff wrt T
947 mutable vector<double> m_dlnActCoeffMolaldT;
948
949 //! 2nd Derivative of log act coeff wrt T
950 mutable vector<double> m_d2lnActCoeffMolaldT2;
951
952 //! Derivative of log act coeff wrt P
953 mutable vector<double> m_dlnActCoeffMolaldP;
954
955private:
956 //! Calculate the log activity coefficients
957 /*!
958 * This function updates the internally stored natural logarithm of the
959 * molality activity coefficients. This is the main routine for
960 * implementing the activity coefficient formulation.
961 */
962 void s_update_lnMolalityActCoeff() const;
963
964 //! Calculation of temperature derivative of activity coefficient
965 /*!
966 * Using internally stored values, this function calculates the temperature
967 * derivative of the logarithm of the activity coefficient for all species
968 * in the mechanism.
969 *
970 * We assume that the activity coefficients are current in this routine. The
971 * solvent activity coefficient is on the molality scale. Its derivative is
972 * too.
973 */
975
976 //! Calculate the temperature 2nd derivative of the activity coefficient
977 /*!
978 * Using internally stored values, this function calculates the temperature
979 * 2nd derivative of the logarithm of the activity coefficient for all
980 * species in the mechanism.
981 *
982 * We assume that the activity coefficients are current in this routine.
983 * Solvent activity coefficient is on the molality scale. Its derivatives
984 * are too.
985 */
987
988 //! Calculate the pressure derivative of the activity coefficient
989 /*!
990 * Using internally stored values, this function calculates the pressure
991 * derivative of the logarithm of the activity coefficient for all species
992 * in the mechanism.
993 *
994 * We assume that the activity coefficients, molalities, and A_Debye are
995 * current. Solvent activity coefficient is on the molality scale. Its
996 * derivatives are too.
997 */
999};
1000
1001}
1002
1003#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:427
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.
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
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.
vector< double > m_tmpV
vector of size m_kk, used as a temporary holding area.
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.
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
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.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
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 gibbs_mole() const override
Molar Gibbs function. Units: J/kmol.
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:562
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:564