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