Cantera 2.6.0
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 * and where we can break down the Gibbs free energy contributions into enthalpy and entropy contributions
55 *
56 * \f[
57 * H^E_i = n X_{Ai} X_{Bi} \sum_m \left( H^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
58 * \f]
59 *
60 * \f[
61 * S^E_i = n X_{Ai} X_{Bi} \sum_m \left( S^{i}_m {\left( X_{Ai} - X_{Bi} \right)}^m \right)
62 * \f]
63 *
64 * where n is the total moles in the solution. The activity of a species defined
65 * in the phase is given by an excess Gibbs free energy formulation.
66 *
67 * \f[
68 * a_k = \gamma_k X_k
69 * \f]
70 *
71 * where
72 *
73 * \f[
74 * R T \ln( \gamma_k )= \frac{d(n G^E)}{d(n_k)}\Bigg|_{n_i}
75 * \f]
76 *
77 * Taking the derivatives results in the following expression
78 * \f[
79 * 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)
80 * + \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)
81 * \f]
82 *
83 * This object inherits from the class VPStandardStateTP. Therefore, the
84 * specification and calculation of all standard state and reference state
85 * values are handled at that level. Various functional forms for the standard
86 * state are permissible. The chemical potential for species *k* is equal to
87 *
88 * \f[
89 * \mu_k(T,P) = \mu^o_k(T, P) + R T \ln(\gamma_k X_k)
90 * \f]
91 *
92 * The partial molar entropy for species *k* is given by the following relation,
93 *
94 * \f[
95 * \tilde{s}_k(T,P) = s^o_k(T,P) - R \ln( \gamma_k X_k )
96 * - R T \frac{d \ln(\gamma_k) }{dT}
97 * \f]
98 *
99 * The partial molar enthalpy for species *k* is given by
100 *
101 * \f[
102 * \tilde{h}_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
103 * \f]
104 *
105 * The partial molar volume for species *k* is
106 *
107 * \f[
108 * \tilde V_k(T,P) = V^o_k(T,P) + R T \frac{d \ln(\gamma_k) }{dP}
109 * \f]
110 *
111 * The partial molar Heat Capacity for species *k* is
112 *
113 * \f[
114 * \tilde{C}_{p,k}(T,P) = C^o_{p,k}(T,P) - 2 R T \frac{d \ln( \gamma_k )}{dT}
115 * - R T^2 \frac{d^2 \ln(\gamma_k) }{{dT}^2}
116 * \f]
117 *
118 * ## %Application within Kinetics Managers
119 *
120 * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / C^s_k, \f$ where
121 * \f$ C^s_k \f$ is a standard concentration defined below and \f$ a_k \f$ are
122 * activities used in the thermodynamic functions. These activity (or
123 * generalized) concentrations are used by kinetics manager classes to compute
124 * the forward and reverse rates of elementary reactions. The activity
125 * concentration,\f$ C^a_k \f$,is given by the following expression.
126 *
127 * \f[
128 * C^a_k = C^s_k X_k = \frac{P}{R T} X_k
129 * \f]
130 *
131 * The standard concentration for species *k* is independent of *k* and equal to
132 *
133 * \f[
134 * C^s_k = C^s = \frac{P}{R T}
135 * \f]
136 *
137 * For example, a bulk-phase binary gas reaction between species j and k,
138 * producing a new gas species l would have the following equation for its rate
139 * of progress variable, \f$ R^1 \f$, which has units of kmol m-3 s-1.
140 *
141 * \f[
142 * R^1 = k^1 C_j^a C_k^a = k^1 (C^s a_j) (C^s a_k)
143 * \f]
144 * where
145 * \f[
146 * C_j^a = C^s a_j \mbox{\quad and \quad} C_k^a = C^s a_k
147 * \f]
148 *
149 * \f$ C_j^a \f$ is the activity concentration of species j, and \f$ C_k^a \f$
150 * is the activity concentration of species k. \f$ C^s \f$ is the standard
151 * concentration. \f$ a_j \f$ is the activity of species j which is equal to the
152 * mole fraction of j.
153 *
154 * The reverse rate constant can then be obtained from the law of microscopic
155 * reversibility and the equilibrium expression for the system.
156 *
157 * \f[
158 * \frac{a_j a_k}{ a_l} = K_a^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} )
159 * \f]
160 *
161 * \f$ K_a^{o,1} \f$ is the dimensionless form of the equilibrium constant,
162 * associated with the pressure dependent standard states \f$ \mu^o_l(T,P) \f$
163 * and their associated activities, \f$ a_l \f$, repeated here:
164 *
165 * \f[
166 * \mu_l(T,P) = \mu^o_l(T, P) + R T \log(a_l)
167 * \f]
168 *
169 * We can switch over to expressing the equilibrium constant in terms of the
170 * reference state chemical potentials
171 *
172 * \f[
173 * K_a^{o,1} = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{P}
174 * \f]
175 *
176 * The concentration equilibrium constant, \f$ K_c \f$, may be obtained by
177 * changing over to activity concentrations. When this is done:
178 *
179 * \f[
180 * \frac{C^a_j C^a_k}{ C^a_l} = C^o K_a^{o,1} = K_c^1 =
181 * \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} ) * \frac{P_{ref}}{RT}
182 * \f]
183 *
184 * %Kinetics managers will calculate the concentration equilibrium constant, \f$
185 * K_c \f$, using the second and third part of the above expression as a
186 * definition for the concentration equilibrium constant.
187 *
188 * For completeness, the pressure equilibrium constant may be obtained as well
189 *
190 * \f[
191 * \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} )
192 * \f]
193 *
194 * \f$ K_p \f$ is the simplest form of the equilibrium constant for ideal gases.
195 * However, it isn't necessarily the simplest form of the equilibrium constant
196 * for other types of phases; \f$ K_c \f$ is used instead because it is
197 * completely general.
198 *
199 * The reverse rate of progress may be written down as
200 * \f[
201 * R^{-1} = k^{-1} C_l^a = k^{-1} (C^o a_l)
202 * \f]
203 *
204 * where we can use the concept of microscopic reversibility to write the
205 * reverse rate constant in terms of the forward rate constant and the
206 * concentration equilibrium constant, \f$ K_c \f$.
207 *
208 * \f[
209 * k^{-1} = k^1 K^1_c
210 * \f]
211 *
212 * \f$k^{-1} \f$ has units of s-1.
213 *
214 * @ingroup thermoprops
215 */
217{
218public:
219 //! Construct a RedlichKisterVPSSTP object from an input file
220 /*!
221 * @param inputFile Name of the input file containing the phase definition.
222 * If blank, an empty phase will be created.
223 * @param id name (ID) of the phase in the input file. If empty, the
224 * first phase definition in the input file will be used.
225 */
226 explicit RedlichKisterVPSSTP(const std::string& inputFile="",
227 const std::string& id="");
228
229 //! Construct and initialize a RedlichKisterVPSSTP ThermoPhase object
230 //! directly from an XML database
231 /*!
232 * @param phaseRef XML phase node containing the description of the phase
233 * @param id id attribute containing the name of the phase.
234 * (default is the empty string)
235 *
236 * @deprecated The XML input format is deprecated and will be removed in
237 * Cantera 3.0.
238 */
239 RedlichKisterVPSSTP(XML_Node& phaseRef, const std::string& id = "");
240
241 virtual std::string type() const {
242 return "RedlichKister";
243 }
244
245 //! @name Molar Thermodynamic Properties
246 //! @{
247
248 virtual doublereal enthalpy_mole() const;
249 virtual doublereal entropy_mole() const;
250 virtual doublereal cp_mole() const;
251 virtual doublereal cv_mole() const;
252
253 /**
254 * @}
255 * @name Activities, Standard States, and Activity Concentrations
256 *
257 * The activity \f$a_k\f$ of a species in solution is
258 * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
259 * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
260 * the chemical potential at unit activity, which depends only
261 * on temperature and pressure.
262 * @{
263 */
264
265 virtual void getLnActivityCoefficients(doublereal* lnac) const;
266
267 //! @}
268 //! @name Partial Molar Properties of the Solution
269 //! @{
270
271 virtual void getChemPotentials(doublereal* mu) const;
272
273 //! Returns an array of partial molar enthalpies for the species in the
274 //! mixture.
275 /*!
276 * Units (J/kmol)
277 *
278 * For this phase, the partial molar enthalpies are equal to the standard
279 * state enthalpies modified by the derivative of the molality-based
280 * activity coefficient wrt temperature
281 *
282 * \f[
283 * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
284 * \f]
285 *
286 * @param hbar Vector of returned partial molar enthalpies
287 * (length m_kk, units = J/kmol)
288 */
289 virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
290
291 //! Returns an array of partial molar entropies for the species in the
292 //! mixture.
293 /*!
294 * Units (J/kmol)
295 *
296 * For this phase, the partial molar enthalpies are equal to the standard
297 * state enthalpies modified by the derivative of the activity coefficient
298 * wrt temperature
299 *
300 * \f[
301 * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
302 * - R \ln( \gamma_k X_k)
303 * - R T \frac{d \ln(\gamma_k) }{dT}
304 * \f]
305 *
306 * @param sbar Vector of returned partial molar entropies
307 * (length m_kk, units = J/kmol/K)
308 */
309 virtual void getPartialMolarEntropies(doublereal* sbar) const;
310
311 //! Returns an array of partial molar entropies for the species in the
312 //! mixture.
313 /*!
314 * Units (J/kmol)
315 *
316 * For this phase, the partial molar enthalpies are equal to the standard
317 * state enthalpies modified by the derivative of the activity coefficient
318 * wrt temperature
319 *
320 * \f[
321 * ???????????????
322 * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
323 * - R \ln( \gamma_k X_k)
324 * - R T \frac{d \ln(\gamma_k) }{dT}
325 * ???????????????
326 * \f]
327 *
328 * @param cpbar Vector of returned partial molar heat capacities
329 * (length m_kk, units = J/kmol/K)
330 */
331 virtual void getPartialMolarCp(doublereal* cpbar) const;
332
333 virtual void getPartialMolarVolumes(doublereal* vbar) const;
334 //! @}
335
336 //! Get the array of temperature second derivatives of the log activity
337 //! coefficients
338 /*!
339 * units = 1/Kelvin
340 *
341 * @param d2lnActCoeffdT2 Output vector of temperature 2nd derivatives of
342 * the log Activity Coefficients. length = m_kk
343 */
344 virtual void getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) const;
345
346 virtual void getdlnActCoeffdT(doublereal* dlnActCoeffdT) const;
347
348 /// @name Initialization
349 /// The following methods are used in the process of constructing
350 /// the phase and setting its parameters from a specification in an
351 /// input file. They are not normally used in application programs.
352 /// To see how they are used, see importPhase().
353
354 virtual void initThermo();
355 virtual void getParameters(AnyMap& phaseNode) const;
356 virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
357
358 //! Add a binary species interaction with the specified parameters
359 /*!
360 * @param speciesA name of the first species
361 * @param speciesB name of the second species
362 * @param excess_enthalpy coefficients of the excess enthalpy polynomial
363 * @param n_enthalpy number of excess enthalpy polynomial coefficients
364 * @param excess_entropy coefficients of the excess entropy polynomial
365 * @param n_entropy number of excess entropy polynomial coefficients
366 */
367 void addBinaryInteraction(const std::string& speciesA, const std::string& speciesB,
368 const double* excess_enthalpy, size_t n_enthalpy,
369 const double* excess_entropy, size_t n_entropy);
370
371 //! @name Derivatives of Thermodynamic Variables needed for Applications
372 //! @{
373
374 virtual void getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds, doublereal* dlnActCoeffds) const;
375 virtual void getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const;
376 virtual void getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const;
377 virtual void getdlnActCoeffdlnN(const size_t ld, doublereal* const dlnActCoeffdlnN);
378 //! @}
379
380private:
381 //! Process an XML node called "binaryNeutralSpeciesParameters"
382 /*!
383 * This node contains all of the parameters necessary to describe the
384 * Redlich-Kister model for a particular binary interaction. This function
385 * reads the XML file and writes the coefficients it finds to an internal
386 * data structures.
387 *
388 * @param xmlBinarySpecies Reference to the XML_Node named
389 * "binaryNeutralSpeciesParameters" containing the binary interaction
390 */
391 void readXMLBinarySpecies(XML_Node& xmlBinarySpecies);
392
393 //! Initialize lengths of local variables after all species have been
394 //! identified.
395 void initLengths();
396
397 //! Update the activity coefficients
398 /*!
399 * This function will be called to update the internally stored natural
400 * logarithm of the activity coefficients
401 */
402 void s_update_lnActCoeff() const;
403
404 //! Update the derivative of the log of the activity coefficients wrt T
405 /*!
406 * This function will be called to update the internally stored derivative
407 * of the natural logarithm of the activity coefficients wrt temperature.
408 */
409 void s_update_dlnActCoeff_dT() const;
410
411 //! Internal routine that calculates the derivative of the activity
412 //! coefficients wrt the mole fractions.
413 /*!
414 * This routine calculates the the derivative of the activity coefficients
415 * wrt to mole fraction with all other mole fractions held constant. This is
416 * strictly not permitted. However, if the resulting matrix is multiplied by
417 * a permissible deltaX vector then everything is ok.
418 *
419 * This is the natural way to handle concentration derivatives in this
420 * routine.
421 */
422 void s_update_dlnActCoeff_dX_() const;
423
424 //! Internal routine that calculates the total derivative of the activity
425 //! coefficients with respect to the log of the mole fractions.
426 /*!
427 * This function will be called to update the internally stored vector of
428 * the total derivatives (that is, not assuming other mole fractions are
429 * constant) of the natural logarithm of the activity coefficients with
430 * respect to the log of the mole fraction.
431 */
433
434protected:
435 //! number of binary interaction expressions
437
438 //! vector of species indices representing species A in the interaction
439 /*!
440 * Each Redlich-Kister excess Gibbs free energy term involves two species,
441 * A and B. This vector identifies species A.
442 */
443 std::vector<size_t> m_pSpecies_A_ij;
444
445 //! vector of species indices representing species B in the interaction
446 /*!
447 * Each Redlich-Kister excess Gibbs free energy term involves two species,
448 * A and B. This vector identifies species B.
449 */
450 std::vector<size_t> m_pSpecies_B_ij;
451
452 //! Vector of the length of the polynomial for the interaction.
453 std::vector<size_t> m_N_ij;
454
455 //! Enthalpy term for the binary mole fraction interaction of the excess
456 //! Gibbs free energy expression
457 mutable std::vector< vector_fp> m_HE_m_ij;
458
459 //! Entropy term for the binary mole fraction interaction of the excess
460 //! Gibbs free energy expression
461 mutable std::vector< vector_fp> m_SE_m_ij;
462
463 //! form of the RedlichKister interaction expression. Currently there is
464 //! only one form.
466
467 //! form of the temperature dependence of the Redlich-Kister interaction
468 //! expression. Currently there is only one form -> constant wrt
469 //! temperature.
471
472 //! Two dimensional array of derivatives of activity coefficients wrt mole
473 //! fractions
475};
476
477}
478
479#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:399
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:30
RedlichKisterVPSSTP is a derived class of GibbsExcessVPSSTP that employs the Redlich-Kister approxima...
void s_update_dlnActCoeff_dX_() const
Internal routine that calculates the derivative of the activity coefficients wrt the mole fractions.
virtual void getParameters(AnyMap &phaseNode) const
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
Array2D dlnActCoeff_dX_
Two dimensional array of derivatives of activity coefficients wrt mole fractions.
virtual void getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const
Get the array of temperature second derivatives of the log activity coefficients.
virtual void getLnActivityCoefficients(doublereal *lnac) const
Get the array of non-dimensional molar-based ln activity coefficients at the current solution tempera...
int formRedlichKister_
form of the RedlichKister interaction expression.
int formTempModel_
form of the temperature dependence of the Redlich-Kister interaction expression.
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies for the species in the mixture.
RedlichKisterVPSSTP(const std::string &inputFile="", const std::string &id="")
Construct a RedlichKisterVPSSTP object from an input file.
virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) const
Get the array of temperature derivatives of the log activity coefficients.
size_t numBinaryInteractions_
number of binary interaction expressions
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
std::vector< vector_fp > m_HE_m_ij
Enthalpy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
virtual void getdlnActCoeffdlnN(const size_t ld, doublereal *const dlnActCoeffdlnN)
Get the array of derivatives of the log activity coefficients with respect to the log of the species ...
void readXMLBinarySpecies(XML_Node &xmlBinarySpecies)
Process an XML node called "binaryNeutralSpeciesParameters".
std::vector< size_t > m_N_ij
Vector of the length of the polynomial for the interaction.
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
virtual void getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const
Get the array of ln mole fraction derivatives of the log activity coefficients - diagonal component o...
virtual void getPartialMolarCp(doublereal *cpbar) const
Returns an array of partial molar entropies for the species in the mixture.
std::vector< size_t > m_pSpecies_A_ij
vector of species indices representing species A in the interaction
void s_update_dlnActCoeff_dT() const
Update the derivative of the log of the activity coefficients wrt T.
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
virtual std::string type() const
String indicating the thermodynamic model implemented.
virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const
Get the array of log species mole number derivatives of the log activity coefficients.
virtual void getdlnActCoeffds(const doublereal dTds, const doublereal *const dXds, doublereal *dlnActCoeffds) const
Get the change in activity coefficients wrt changes in state (temp, mole fraction,...
std::vector< size_t > m_pSpecies_B_ij
vector of species indices representing species B in the interaction
void initLengths()
Initialize lengths of local variables after all species have been identified.
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.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
std::vector< vector_fp > m_SE_m_ij
Entropy term for the binary mole fraction interaction of the excess Gibbs free energy expression.
void addBinaryInteraction(const std::string &speciesA, const std::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.
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:103
Namespace for the Cantera kernel.
Definition: AnyMap.h:29