Loading [MathJax]/extensions/tex2jax.js
Cantera  3.2.0a1
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
RedlichKisterVPSSTP.h
Go to the documentation of this file.
1/**
2 * @file RedlichKisterVPSSTP.h (see @ref thermoprops and class @link
3 * Cantera::RedlichKisterVPSSTP RedlichKisterVPSSTP@endlink).
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
9#ifndef CT_REDLICHKISTERVPSSTP_H
10#define CT_REDLICHKISTERVPSSTP_H
11
13
14namespace Cantera
15{
16
17//! RedlichKisterVPSSTP is a derived class of GibbsExcessVPSSTP that employs the
18//! Redlich-Kister approximation for the excess Gibbs free energy
19/*!
20 * RedlichKisterVPSSTP derives from class GibbsExcessVPSSTP which is derived
21 * from VPStandardStateTP, and overloads the virtual methods defined there with
22 * ones that use expressions appropriate for the Redlich Kister Excess Gibbs
23 * free energy approximation.
24 *
25 * The independent unknowns are pressure, temperature, and mass fraction.
26 *
27 * ## Specification of Species Standard State Properties
28 *
29 * All species are defined to have standard states that depend upon both the
30 * temperature and the pressure. The Redlich-Kister approximation assumes
31 * symmetric standard states, where all of the standard state assume that the
32 * species are in pure component states at the temperature and pressure of the
33 * solution. I don't think it prevents, however, some species from being dilute
34 * in the solution.
35 *
36 * ## Specification of Solution Thermodynamic Properties
37 *
38 * The molar excess Gibbs free energy is given by the following formula which is
39 * a sum over interactions *i*. Each of the interactions are binary interactions
40 * involving two of the species in the phase, denoted, *Ai* and *Bi*. This is
41 * the generalization of the Redlich-Kister formulation for a phase that has
42 * more than 2 species.
43 *
44 * @f[
45 * G^E = \sum_{i} G^E_{i}
46 * @f]
47 *
48 * where
49 *
50 * @f[
51 * G^E_{i} = n X_{Ai} X_{Bi} \sum_m \left( A^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
52 * @f]
53 *
54 * where n is the total moles in the solution and where we can break down the Gibbs free
55 * energy contributions into enthalpy and entropy contributions by defining
56 * @f$ A^i_m = H^i_m - T S^i_m @f$ :
57 *
58 * @f[
59 * H^E_i = n X_{Ai} X_{Bi} \sum_m \left( H^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
60 * @f]
61 *
62 * @f[
63 * S^E_i = n X_{Ai} X_{Bi} \sum_m \left( S^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
64 * @f]
65 *
66 * The activity of a species defined in the phase is given by an excess Gibbs free
67 * energy formulation:
68 *
69 * @f[
70 * a_k = \gamma_k X_k
71 * @f]
72 *
73 * where
74 *
75 * @f[
76 * R T \ln( \gamma_k )= \frac{d(n G^E)}{d(n_k)}\Bigg|_{n_i}
77 * @f]
78 *
79 * Taking the derivatives results in the following expression
80 * @f[
81 * R T \ln( \gamma_k )= \sum_i \delta_{Ai,k} (1 - X_{Ai}) X_{Bi} \sum_m \left( A^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
82 * + \sum_i \delta_{Ai,k} X_{Ai} X_{Bi} \sum_m \left( A^{i}_0 + A^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^{m-1} (1 - X_{Ai} + X_{Bi}) \right)
83 * @f]
84 *
85 * Evaluating thermodynamic properties requires the following derivatives of
86 * @f$ \ln(\gamma_k) @f$:
87 *
88 * @f[
89 * \frac{d \ln( \gamma_k )}{dT} = - \frac{1}{RT^2} \left( \sum_i \delta_{Ai,k} (1 - X_{Ai}) X_{Bi} \sum_m \left( H^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
90 * + \sum_i \delta_{Ai,k} X_{Ai} X_{Bi} \sum_m \left( H^{i}_0 + H^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^{m-1} (1 - X_{Ai} + X_{Bi}) \right) \right)
91 * @f]
92 *
93 * and
94 *
95 * @f[
96 * \frac{d^2 \ln( \gamma_k )}{dT^2} = -\frac{2}{T} \frac{d \ln( \gamma_k )}{dT}
97 * @f]
98 *
99 * This object inherits from the class VPStandardStateTP. Therefore, the
100 * specification and calculation of all standard state and reference state
101 * values are handled at that level. Various functional forms for the standard
102 * state are permissible. The chemical potential for species *k* is equal to
103 *
104 * @f[
105 * \mu_k(T,P) = \mu^o_k(T, P) + R T \ln(\gamma_k X_k)
106 * @f]
107 *
108 * The partial molar entropy for species *k* is given by the following relation,
109 *
110 * @f[
111 * \tilde{s}_k(T,P) = s^o_k(T,P) - R \ln( \gamma_k X_k )
112 * - R T \frac{d \ln(\gamma_k) }{dT}
113 * @f]
114 *
115 * The partial molar enthalpy for species *k* is given by
116 *
117 * @f[
118 * \tilde{h}_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
119 * @f]
120 *
121 * The partial molar volume for species *k* is
122 *
123 * @f[
124 * \tilde V_k(T,P) = V^o_k(T,P) + R T \frac{d \ln(\gamma_k) }{dP}
125 * @f]
126 *
127 * The partial molar Heat Capacity for species *k* is
128 *
129 * @f[
130 * \tilde{C}_{p,k}(T,P) = C^o_{p,k}(T,P) - 2 R T \frac{d \ln( \gamma_k )}{dT}
131 * - R T^2 \frac{d^2 \ln(\gamma_k) }{{dT}^2} = C^o_{p,k}(T,P)
132 * @f]
133 *
134 * ## Application within Kinetics Managers
135 *
136 * @f$ C^a_k @f$ are defined such that @f$ a_k = C^a_k / C^s_k, @f$ where
137 * @f$ C^s_k @f$ is a standard concentration defined below and @f$ a_k @f$ are
138 * activities used in the thermodynamic functions. These activity (or
139 * generalized) concentrations are used by kinetics manager classes to compute
140 * the forward and reverse rates of elementary reactions. The activity
141 * concentration,@f$ C^a_k @f$,is given by the following expression.
142 *
143 * @f[
144 * C^a_k = C^s_k X_k = \frac{P}{R T} X_k
145 * @f]
146 *
147 * The standard concentration for species *k* is independent of *k* and equal to
148 *
149 * @f[
150 * C^s_k = C^s = \frac{P}{R T}
151 * @f]
152 *
153 * For example, a bulk-phase binary gas reaction between species j and k,
154 * producing a new gas species l would have the following equation for its rate
155 * of progress variable, @f$ R^1 @f$, which has units of kmol m-3 s-1.
156 *
157 * @f[
158 * R^1 = k^1 C_j^a C_k^a = k^1 (C^s a_j) (C^s a_k)
159 * @f]
160 * where
161 * @f[
162 * C_j^a = C^s a_j \mbox{\quad and \quad} C_k^a = C^s a_k
163 * @f]
164 *
165 * @f$ C_j^a @f$ is the activity concentration of species j, and @f$ C_k^a @f$
166 * is the activity concentration of species k. @f$ C^s @f$ is the standard
167 * concentration. @f$ a_j @f$ is the activity of species j which is equal to the
168 * mole fraction of j.
169 *
170 * The reverse rate constant can then be obtained from the law of microscopic
171 * reversibility and the equilibrium expression for the system.
172 *
173 * @f[
174 * \frac{a_j a_k}{ a_l} = K_a^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} )
175 * @f]
176 *
177 * @f$ K_a^{o,1} @f$ is the dimensionless form of the equilibrium constant,
178 * associated with the pressure dependent standard states @f$ \mu^o_l(T,P) @f$
179 * and their associated activities, @f$ a_l @f$, repeated here:
180 *
181 * @f[
182 * \mu_l(T,P) = \mu^o_l(T, P) + R T \ln a_l
183 * @f]
184 *
185 * We can switch over to expressing the equilibrium constant in terms of the
186 * reference state chemical potentials
187 *
188 * @f[
189 * K_a^{o,1} = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{P}
190 * @f]
191 *
192 * The concentration equilibrium constant, @f$ K_c @f$, may be obtained by
193 * changing over to activity concentrations. When this is done:
194 *
195 * @f[
196 * \frac{C^a_j C^a_k}{ C^a_l} = C^o K_a^{o,1} = K_c^1 =
197 * \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{RT}
198 * @f]
199 *
200 * %Kinetics managers will calculate the concentration equilibrium constant, @f$
201 * K_c @f$, using the second and third part of the above expression as a
202 * definition for the concentration equilibrium constant.
203 *
204 * For completeness, the pressure equilibrium constant may be obtained as well
205 *
206 * @f[
207 * \frac{P_j P_k}{ P_l P_{ref}} = K_p^1 = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} )
208 * @f]
209 *
210 * @f$ K_p @f$ is the simplest form of the equilibrium constant for ideal gases.
211 * However, it isn't necessarily the simplest form of the equilibrium constant
212 * for other types of phases; @f$ K_c @f$ is used instead because it is
213 * completely general.
214 *
215 * The reverse rate of progress may be written down as
216 * @f[
217 * R^{-1} = k^{-1} C_l^a = k^{-1} (C^o a_l)
218 * @f]
219 *
220 * where we can use the concept of microscopic reversibility to write the
221 * reverse rate constant in terms of the forward rate constant and the
222 * concentration equilibrium constant, @f$ K_c @f$.
223 *
224 * @f[
225 * k^{-1} = k^1 K^1_c
226 * @f]
227 *
228 * @f$ k^{-1} @f$ has units of s-1.
229 *
230 * @ingroup thermoprops
231 */
233{
234public:
235 //! Construct a RedlichKisterVPSSTP object from an input file
236 /*!
237 * @param inputFile Name of the input file containing the phase definition.
238 * If blank, an empty phase will be created.
239 * @param id name (ID) of the phase in the input file. If empty, the
240 * first phase definition in the input file will be used.
241 */
242 explicit RedlichKisterVPSSTP(const string& inputFile="", const string& id="");
243
244 string type() const override {
245 return "Redlich-Kister";
246 }
247
248 //! @name Molar Thermodynamic Properties
249 //! @{
250
251 double cv_mole() const override;
252
253 //! @}
254 //! @name Activities, Standard States, and Activity Concentrations
255 //!
256 //! The activity @f$ a_k @f$ of a species in solution is
257 //! related to the chemical potential by @f[ \mu_k = \mu_k^0(T)
258 //! + \hat R T \ln a_k. @f] The quantity @f$ \mu_k^0(T,P) @f$ is
259 //! the chemical potential at unit activity, which depends only
260 //! on temperature and pressure.
261 //! @{
262
263 void getLnActivityCoefficients(double* lnac) const override;
264
265 //! @}
266 //! @name Partial Molar Properties of the Solution
267 //! @{
268
269 void getChemPotentials(double* mu) const override;
270
271 //! Returns an array of partial molar enthalpies for the species in the
272 //! mixture.
273 /*!
274 * Units (J/kmol)
275 *
276 * For this phase, the partial molar enthalpies are equal to the standard
277 * state enthalpies modified by the derivative of the molality-based
278 * activity coefficient wrt temperature
279 *
280 * @f[
281 * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
282 * @f]
283 *
284 * @param hbar Vector of returned partial molar enthalpies
285 * (length m_kk, units = J/kmol)
286 */
287 void getPartialMolarEnthalpies(double* hbar) const override;
288
289 //! Returns an array of partial molar entropies for the species in the
290 //! mixture.
291 /*!
292 * For this phase, the partial molar entropies are equal to the standard
293 * state entropies modified by the derivative of the activity coefficient
294 * with respect to temperature:
295 *
296 * @f[
297 * \bar s_k(T,P) = s^o_k(T,P) - R \ln( \gamma_k X_k)
298 * - R T \frac{d \ln(\gamma_k) }{dT}
299 * @f]
300 *
301 * @param sbar Vector of returned partial molar entropies
302 * (length m_kk, units = J/kmol/K)
303 */
304 void getPartialMolarEntropies(double* sbar) const override;
305
306 //! Returns an array of partial molar heat capacities for the species in the
307 //! mixture.
308 /*!
309 * Units (J/kmol/K)
310 *
311 * For this phase, the partial molar heat capacities are equal to the standard
312 * state heat capacities:
313 *
314 * @f[
315 * \tilde{C}_{p,k}(T,P) = C^o_{p,k}(T,P)
316 * @f]
317 *
318 * @param cpbar Vector of returned partial molar heat capacities
319 * (length m_kk, units = J/kmol/K)
320 */
321 void getPartialMolarCp(double* cpbar) const override;
322
323 void getPartialMolarVolumes(double* vbar) const override;
324 //! @}
325
326 //! Get the array of temperature second derivatives of the log activity
327 //! coefficients
328 /*!
329 * units = 1/Kelvin
330 *
331 * @param d2lnActCoeffdT2 Output vector of temperature 2nd derivatives of
332 * the log Activity Coefficients. length = m_kk
333 */
334 void getd2lnActCoeffdT2(double* d2lnActCoeffdT2) const;
335
336 void getdlnActCoeffdT(double* dlnActCoeffdT) const override;
337
338 //! @name Initialization
339 //!
340 //! The following methods are used in the process of constructing
341 //! the phase and setting its parameters from a specification in an
342 //! input file. They are not normally used in application programs.
343 //! To see how they are used, see importPhase().
344
345 void initThermo() override;
346 void getParameters(AnyMap& phaseNode) const override;
347
348 //! Add a binary species interaction with the specified parameters
349 /*!
350 * @param speciesA name of the first species
351 * @param speciesB name of the second species
352 * @param excess_enthalpy coefficients of the excess enthalpy polynomial
353 * @param n_enthalpy number of excess enthalpy polynomial coefficients
354 * @param excess_entropy coefficients of the excess entropy polynomial
355 * @param n_entropy number of excess entropy polynomial coefficients
356 */
357 void addBinaryInteraction(const string& speciesA, const string& speciesB,
358 const double* excess_enthalpy, size_t n_enthalpy,
359 const double* excess_entropy, size_t n_entropy);
360
361 //! @name Derivatives of Thermodynamic Variables needed for Applications
362 //! @{
363
364 void getdlnActCoeffds(const double dTds, const double* const dXds,
365 double* dlnActCoeffds) const override;
366 void getdlnActCoeffdlnX_diag(double* dlnActCoeffdlnX_diag) const override;
367 void getdlnActCoeffdlnN_diag(double* dlnActCoeffdlnN_diag) const override;
368 void getdlnActCoeffdlnN(const size_t ld, double* const dlnActCoeffdlnN) override;
369 //! @}
370
371private:
372 //! Initialize lengths of local variables after all species have been
373 //! identified.
374 void initLengths();
375
376 //! Update the activity coefficients
377 /*!
378 * This function will be called to update the internally stored natural
379 * logarithm of the activity coefficients
380 */
381 void s_update_lnActCoeff() const;
382
383 //! Update the derivative of the log of the activity coefficients wrt T
384 /*!
385 * This function will be called to update the internally stored derivative
386 * of the natural logarithm of the activity coefficients wrt temperature.
387 */
388 void s_update_dlnActCoeff_dT() const;
389
390 //! Internal routine that calculates the derivative of the activity
391 //! coefficients wrt the mole fractions.
392 /*!
393 * This routine calculates the the derivative of the activity coefficients
394 * wrt to mole fraction with all other mole fractions held constant. This is
395 * strictly not permitted. However, if the resulting matrix is multiplied by
396 * a permissible deltaX vector then everything is ok.
397 *
398 * This is the natural way to handle concentration derivatives in this
399 * routine.
400 */
401 void s_update_dlnActCoeff_dX_() const;
402
403 //! Internal routine that calculates the total derivative of the activity
404 //! coefficients with respect to the log of the mole fractions.
405 /*!
406 * This function will be called to update the internally stored vector of
407 * the total derivatives (that is, not assuming other mole fractions are
408 * constant) of the natural logarithm of the activity coefficients with
409 * respect to the log of the mole fraction.
410 */
412
413protected:
414 //! vector of species indices representing species A in the interaction
415 /*!
416 * Each Redlich-Kister excess Gibbs free energy term involves two species,
417 * A and B. This vector identifies species A.
418 */
419 vector<size_t> m_pSpecies_A_ij;
420
421 //! vector of species indices representing species B in the interaction
422 /*!
423 * Each Redlich-Kister excess Gibbs free energy term involves two species,
424 * A and B. This vector identifies species B.
425 */
426 vector<size_t> m_pSpecies_B_ij;
427
428 //! Enthalpy term for the binary mole fraction interaction of the excess
429 //! Gibbs free energy expression
430 vector<vector<double>> m_HE_m_ij;
431
432 //! Entropy term for the binary mole fraction interaction of the excess
433 //! Gibbs free energy expression
434 vector<vector<double>> m_SE_m_ij;
435
436 //! Two dimensional array of derivatives of activity coefficients wrt mole
437 //! fractions
439};
440
441}
442
443#endif
Header for intermediate ThermoPhase object for phases which employ Gibbs excess free energy based for...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:432
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
GibbsExcessVPSSTP is a derived class of ThermoPhase that handles variable pressure standard state met...
RedlichKisterVPSSTP is a derived class of GibbsExcessVPSSTP that employs the Redlich-Kister approxima...
void getdlnActCoeffds(const double dTds, const double *const dXds, double *dlnActCoeffds) const override
Get the change in activity coefficients wrt changes in state (temp, mole fraction,...
void s_update_dlnActCoeff_dX_() const
Internal routine that calculates the derivative of the activity coefficients wrt the mole fractions.
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
Array2D dlnActCoeff_dX_
Two dimensional array of derivatives of activity coefficients wrt mole fractions.
vector< vector< double > > m_SE_m_ij
Entropy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
string type() const override
String indicating the thermodynamic model implemented.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
vector< size_t > m_pSpecies_A_ij
vector of species indices representing species A in the interaction
double cv_mole() const override
Molar heat capacity at constant volume. Units: J/kmol/K.
void s_update_dlnActCoeff_dT() const
Update the derivative of the log of the activity coefficients wrt T.
vector< size_t > m_pSpecies_B_ij
vector of species indices representing species B in the interaction
void addBinaryInteraction(const string &speciesA, const string &speciesB, const double *excess_enthalpy, size_t n_enthalpy, const double *excess_entropy, size_t n_entropy)
Add a binary species interaction with the specified parameters.
vector< vector< double > > m_HE_m_ij
Enthalpy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
void getdlnActCoeffdT(double *dlnActCoeffdT) const override
Get the array of temperature derivatives of the log activity coefficients.
void getPartialMolarCp(double *cpbar) const override
Returns an array of partial molar heat capacities for the species in the mixture.
void initLengths()
Initialize lengths of local variables after all species have been identified.
void getLnActivityCoefficients(double *lnac) const override
Get the array of non-dimensional molar-based ln activity coefficients at the current solution tempera...
void s_update_dlnActCoeff_dlnX_diag() const
Internal routine that calculates the total derivative of the activity coefficients with respect to th...
void s_update_lnActCoeff() const
Update the activity coefficients.
void getdlnActCoeffdlnX_diag(double *dlnActCoeffdlnX_diag) const override
Get the array of ln mole fraction derivatives of the log activity coefficients - diagonal component o...
void getdlnActCoeffdlnN_diag(double *dlnActCoeffdlnN_diag) const override
Get the array of log species mole number derivatives of the log activity coefficients.
void getd2lnActCoeffdT2(double *d2lnActCoeffdT2) const
Get the array of temperature second derivatives of the log activity coefficients.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies for the species in the mixture.
void getdlnActCoeffdlnN(const size_t ld, double *const dlnActCoeffdlnN) override
Get the array of derivatives of the log activity coefficients with respect to the log of the species ...
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595