Cantera  3.0.0
Loading...
Searching...
No Matches
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 enthalpy_mole() const override;
252 double entropy_mole() const override;
253 double cp_mole() const override;
254 double cv_mole() const override;
255
256 //! @}
257 //! @name Activities, Standard States, and Activity Concentrations
258 //!
259 //! The activity @f$ a_k @f$ of a species in solution is
260 //! related to the chemical potential by @f[ \mu_k = \mu_k^0(T)
261 //! + \hat R T \ln a_k. @f] The quantity @f$ \mu_k^0(T,P) @f$ is
262 //! the chemical potential at unit activity, which depends only
263 //! on temperature and pressure.
264 //! @{
265
266 void getLnActivityCoefficients(double* lnac) const override;
267
268 //! @}
269 //! @name Partial Molar Properties of the Solution
270 //! @{
271
272 void getChemPotentials(double* mu) const override;
273
274 //! Returns an array of partial molar enthalpies for the species in the
275 //! mixture.
276 /*!
277 * Units (J/kmol)
278 *
279 * For this phase, the partial molar enthalpies are equal to the standard
280 * state enthalpies modified by the derivative of the molality-based
281 * activity coefficient wrt temperature
282 *
283 * @f[
284 * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
285 * @f]
286 *
287 * @param hbar Vector of returned partial molar enthalpies
288 * (length m_kk, units = J/kmol)
289 */
290 void getPartialMolarEnthalpies(double* hbar) const override;
291
292 //! Returns an array of partial molar entropies for the species in the
293 //! mixture.
294 /*!
295 * For this phase, the partial molar entropies are equal to the standard
296 * state entropies modified by the derivative of the activity coefficient
297 * with respect to temperature:
298 *
299 * @f[
300 * \bar s_k(T,P) = s^o_k(T,P) - R \ln( \gamma_k X_k)
301 * - R T \frac{d \ln(\gamma_k) }{dT}
302 * @f]
303 *
304 * @param sbar Vector of returned partial molar entropies
305 * (length m_kk, units = J/kmol/K)
306 */
307 void getPartialMolarEntropies(double* sbar) const override;
308
309 //! Returns an array of partial molar heat capacities for the species in the
310 //! mixture.
311 /*!
312 * Units (J/kmol/K)
313 *
314 * For this phase, the partial molar heat capacities are equal to the standard
315 * state heat capacities:
316 *
317 * @f[
318 * \tilde{C}_{p,k}(T,P) = C^o_{p,k}(T,P)
319 * @f]
320 *
321 * @param cpbar Vector of returned partial molar heat capacities
322 * (length m_kk, units = J/kmol/K)
323 */
324 void getPartialMolarCp(double* cpbar) const override;
325
326 void getPartialMolarVolumes(double* vbar) const override;
327 //! @}
328
329 //! Get the array of temperature second derivatives of the log activity
330 //! coefficients
331 /*!
332 * units = 1/Kelvin
333 *
334 * @param d2lnActCoeffdT2 Output vector of temperature 2nd derivatives of
335 * the log Activity Coefficients. length = m_kk
336 */
337 void getd2lnActCoeffdT2(double* d2lnActCoeffdT2) const;
338
339 void getdlnActCoeffdT(double* dlnActCoeffdT) const override;
340
341 //! @name Initialization
342 //!
343 //! The following methods are used in the process of constructing
344 //! the phase and setting its parameters from a specification in an
345 //! input file. They are not normally used in application programs.
346 //! To see how they are used, see importPhase().
347
348 void initThermo() override;
349 void getParameters(AnyMap& phaseNode) const override;
350
351 //! Add a binary species interaction with the specified parameters
352 /*!
353 * @param speciesA name of the first species
354 * @param speciesB name of the second species
355 * @param excess_enthalpy coefficients of the excess enthalpy polynomial
356 * @param n_enthalpy number of excess enthalpy polynomial coefficients
357 * @param excess_entropy coefficients of the excess entropy polynomial
358 * @param n_entropy number of excess entropy polynomial coefficients
359 */
360 void addBinaryInteraction(const string& speciesA, const string& speciesB,
361 const double* excess_enthalpy, size_t n_enthalpy,
362 const double* excess_entropy, size_t n_entropy);
363
364 //! @name Derivatives of Thermodynamic Variables needed for Applications
365 //! @{
366
367 void getdlnActCoeffds(const double dTds, const double* const dXds,
368 double* dlnActCoeffds) const override;
369 void getdlnActCoeffdlnX_diag(double* dlnActCoeffdlnX_diag) const override;
370 void getdlnActCoeffdlnN_diag(double* dlnActCoeffdlnN_diag) const override;
371 void getdlnActCoeffdlnN(const size_t ld, double* const dlnActCoeffdlnN) override;
372 //! @}
373
374private:
375 //! Initialize lengths of local variables after all species have been
376 //! identified.
377 void initLengths();
378
379 //! Update the activity coefficients
380 /*!
381 * This function will be called to update the internally stored natural
382 * logarithm of the activity coefficients
383 */
384 void s_update_lnActCoeff() const;
385
386 //! Update the derivative of the log of the activity coefficients wrt T
387 /*!
388 * This function will be called to update the internally stored derivative
389 * of the natural logarithm of the activity coefficients wrt temperature.
390 */
391 void s_update_dlnActCoeff_dT() const;
392
393 //! Internal routine that calculates the derivative of the activity
394 //! coefficients wrt the mole fractions.
395 /*!
396 * This routine calculates the the derivative of the activity coefficients
397 * wrt to mole fraction with all other mole fractions held constant. This is
398 * strictly not permitted. However, if the resulting matrix is multiplied by
399 * a permissible deltaX vector then everything is ok.
400 *
401 * This is the natural way to handle concentration derivatives in this
402 * routine.
403 */
404 void s_update_dlnActCoeff_dX_() const;
405
406 //! Internal routine that calculates the total derivative of the activity
407 //! coefficients with respect to the log of the mole fractions.
408 /*!
409 * This function will be called to update the internally stored vector of
410 * the total derivatives (that is, not assuming other mole fractions are
411 * constant) of the natural logarithm of the activity coefficients with
412 * respect to the log of the mole fraction.
413 */
415
416protected:
417 //! vector of species indices representing species A in the interaction
418 /*!
419 * Each Redlich-Kister excess Gibbs free energy term involves two species,
420 * A and B. This vector identifies species A.
421 */
422 vector<size_t> m_pSpecies_A_ij;
423
424 //! vector of species indices representing species B in the interaction
425 /*!
426 * Each Redlich-Kister excess Gibbs free energy term involves two species,
427 * A and B. This vector identifies species B.
428 */
429 vector<size_t> m_pSpecies_B_ij;
430
431 //! Enthalpy term for the binary mole fraction interaction of the excess
432 //! Gibbs free energy expression
433 vector<vector<double>> m_HE_m_ij;
434
435 //! Entropy term for the binary mole fraction interaction of the excess
436 //! Gibbs free energy expression
437 vector<vector<double>> m_SE_m_ij;
438
439 //! Two dimensional array of derivatives of activity coefficients wrt mole
440 //! fractions
442};
443
444}
445
446#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:427
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.
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
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.
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
void getdlnActCoeffdT(double *dlnActCoeffdT) const override
Get the array of temperature derivatives of the log activity coefficients.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
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:564