Cantera  2.3.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 http://www.cantera.org/license.txt for license and copyright information.
8 
9 #ifndef CT_REDLICHKISTERVPSSTP_H
10 #define CT_REDLICHKISTERVPSSTP_H
11 
13 
14 namespace 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 reate 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 {
218 public:
219  //! Constructor
220  /*!
221  * This doesn't do much more than initialize constants with default values.
222  */
224 
225  //! Construct and initialize a RedlichKisterVPSSTP ThermoPhase object
226  //! directly from an XML input file
227  /*!
228  * @param inputFile Name of the input file containing the phase XML data
229  * to set up the object
230  * @param id ID of the phase in the input file. Defaults to the
231  * empty string.
232  */
233  RedlichKisterVPSSTP(const std::string& inputFile, const std::string& id = "");
234 
235  //! Construct and initialize a RedlichKisterVPSSTP ThermoPhase object
236  //! directly from an XML database
237  /*!
238  * @param phaseRef XML phase node containing the description of the phase
239  * @param id id attribute containing the name of the phase.
240  * (default is the empty string)
241  */
242  RedlichKisterVPSSTP(XML_Node& phaseRef, const std::string& id = "");
243 
245  RedlichKisterVPSSTP& operator=(const RedlichKisterVPSSTP& b);
246  virtual ThermoPhase* duplMyselfAsThermoPhase() const;
247 
248  virtual std::string type() const {
249  return "RedlichKister";
250  }
251 
252  //! @name Molar Thermodynamic Properties
253  //! @{
254 
255  virtual doublereal enthalpy_mole() const;
256  virtual doublereal entropy_mole() const;
257  virtual doublereal cp_mole() const;
258  virtual doublereal cv_mole() const;
259 
260  /**
261  * @}
262  * @name Activities, Standard States, and Activity Concentrations
263  *
264  * The activity \f$a_k\f$ of a species in solution is
265  * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
266  * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
267  * the chemical potential at unit activity, which depends only
268  * on temperature and pressure.
269  * @{
270  */
271 
272  virtual void getLnActivityCoefficients(doublereal* lnac) const;
273 
274  //@}
275  /// @name Partial Molar Properties of the Solution
276  //@{
277 
278  virtual void getChemPotentials(doublereal* mu) const;
279 
280  //! Returns an array of partial molar enthalpies for the species in the
281  //! mixture.
282  /*!
283  * Units (J/kmol)
284  *
285  * For this phase, the partial molar enthalpies are equal to the standard
286  * state enthalpies modified by the derivative of the molality-based
287  * activity coefficient wrt temperature
288  *
289  * \f[
290  * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
291  * \f]
292  *
293  * @param hbar Vector of returned partial molar enthalpies
294  * (length m_kk, units = J/kmol)
295  */
296  virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
297 
298  //! Returns an array of partial molar entropies for the species in the
299  //! mixture.
300  /*!
301  * Units (J/kmol)
302  *
303  * For this phase, the partial molar enthalpies are equal to the standard
304  * state enthalpies modified by the derivative of the activity coefficient
305  * wrt temperature
306  *
307  * \f[
308  * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
309  * - R \ln( \gamma_k X_k)
310  * - R T \frac{d \ln(\gamma_k) }{dT}
311  * \f]
312  *
313  * @param sbar Vector of returned partial molar entropies
314  * (length m_kk, units = J/kmol/K)
315  */
316  virtual void getPartialMolarEntropies(doublereal* sbar) const;
317 
318  //! Returns an array of partial molar entropies for the species in the
319  //! mixture.
320  /*!
321  * Units (J/kmol)
322  *
323  * For this phase, the partial molar enthalpies are equal to the standard
324  * state enthalpies modified by the derivative of the activity coefficient
325  * wrt temperature
326  *
327  * \f[
328  * ???????????????
329  * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
330  * - R \ln( \gamma_k X_k)
331  * - R T \frac{d \ln(\gamma_k) }{dT}
332  * ???????????????
333  * \f]
334  *
335  * @param cpbar Vector of returned partial molar heat capacities
336  * (length m_kk, units = J/kmol/K)
337  */
338  virtual void getPartialMolarCp(doublereal* cpbar) const;
339 
340  virtual void getPartialMolarVolumes(doublereal* vbar) const;
341 
342  //! Get the array of temperature second derivatives of the log activity
343  //! coefficients
344  /*!
345  * units = 1/Kelvin
346  *
347  * @param d2lnActCoeffdT2 Output vector of temperature 2nd derivatives of
348  * the log Activity Coefficients. length = m_kk
349  */
350  virtual void getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) const;
351 
352  virtual void getdlnActCoeffdT(doublereal* dlnActCoeffdT) const;
353 
354  /// @}
355  /// @name Initialization
356  /// The following methods are used in the process of constructing
357  /// the phase and setting its parameters from a specification in an
358  /// input file. They are not normally used in application programs.
359  /// To see how they are used, see importPhase().
360 
361  virtual void initThermo();
362  virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
363 
364  //! @}
365  //! @name Derivatives of Thermodynamic Variables needed for Applications
366  //! @{
367 
368  virtual void getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds, doublereal* dlnActCoeffds) const;
369  virtual void getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const;
370  virtual void getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const;
371  virtual void getdlnActCoeffdlnN(const size_t ld, doublereal* const dlnActCoeffdlnN);
372 
373  //@}
374 
375 private:
376  //! Process an XML node called "binaryNeutralSpeciesParameters"
377  /*!
378  * This node contains all of the parameters necessary to describe the
379  * Redlich-Kister model for a particular binary interaction. This function
380  * reads the XML file and writes the coefficients it finds to an internal
381  * data structures.
382  *
383  * @param xmlBinarySpecies Reference to the XML_Node named
384  * "binaryNeutralSpeciesParameters" containing the binary interaction
385  */
386  void readXMLBinarySpecies(XML_Node& xmlBinarySpecies);
387 
388  //! Resize internal arrays within the object that depend upon the number
389  //! of binary Redlich-Kister interaction terms
390  /*!
391  * @param num Number of binary Redlich-Kister interaction terms
392  */
393  void resizeNumInteractions(const size_t num);
394 
395  //! Initialize lengths of local variables after all species have been
396  //! identified.
397  void initLengths();
398 
399  //! Update the activity coefficients
400  /*!
401  * This function will be called to update the internally stored natural
402  * logarithm of the activity coefficients
403  */
404  void s_update_lnActCoeff() const;
405 
406  //! Update the derivative of the log of the activity coefficients wrt T
407  /*!
408  * This function will be called to update the internally stored derivative
409  * of the natural logarithm of the activity coefficients wrt temperature.
410  */
411  void s_update_dlnActCoeff_dT() const;
412 
413  //! Internal routine that calculates the derivative of the activity
414  //! coefficients wrt the mole fractions.
415  /*!
416  * This routine calculates the the derivative of the activity coefficients
417  * wrt to mole fraction with all other mole fractions held constant. This is
418  * strictly not permitted. However, if the resulting matrix is multiplied by
419  * a permissible deltaX vector then everything is ok.
420  *
421  * This is the natural way to handle concentration derivatives in this
422  * routine.
423  */
424  void s_update_dlnActCoeff_dX_() const;
425 
426  //! Internal routine that calculates the total derivative of the activity
427  //! coefficients with respect to the log of the mole fractions.
428  /*!
429  * This function will be called to update the internally stored vector of
430  * the total derivatives (i.e. not assuming other mole fractions are
431  * constant) of the natural logarithm of the activity coefficients with
432  * respect to the log of the mole fraction.
433  */
434  void s_update_dlnActCoeff_dlnX_diag() const;
435 
436 public:
437  //! Utility routine that calculates a literature expression
438  /*!
439  * @param VintOut Output contribution to the voltage corresponding to
440  * nonideal term
441  * @param voltsOut Output contribution to the voltage corresponding to
442  * nonideal term and mf term
443  * @deprecated Probably broken. To be removed after Cantera 2.3.
444  */
445  void Vint(double& VintOut, double& voltsOut);
446 
447 protected:
448  //! number of binary interaction expressions
450 
451  //! vector of species indices representing species A in the interaction
452  /*!
453  * Each Redlich-Kister excess Gibbs free energy term involves two species,
454  * A and B. This vector identifies species A.
455  */
456  std::vector<size_t> m_pSpecies_A_ij;
457 
458  //! vector of species indices representing species B in the interaction
459  /*!
460  * Each Redlich-Kister excess Gibbs free energy term involves two species,
461  * A and B. This vector identifies species B.
462  */
463  std::vector<size_t> m_pSpecies_B_ij;
464 
465  //! Vector of the length of the polynomial for the interaction.
466  std::vector<size_t> m_N_ij;
467 
468  //! Enthalpy term for the binary mole fraction interaction of the excess
469  //! Gibbs free energy expression
470  mutable std::vector< vector_fp> m_HE_m_ij;
471 
472  //! Entropy term for the binary mole fraction interaction of the excess
473  //! Gibbs free energy expression
474  mutable std::vector< vector_fp> m_SE_m_ij;
475 
476  //! form of the RedlichKister interaction expression. Currently there is
477  //! only one form.
479 
480  //! form of the temperature dependence of the Redlich-Kister interaction
481  //! expression. Currently there is only one form -> constant wrt
482  //! temperature.
484 
485  //! Two dimensional array of derivatives of activity coefficients wrt mole
486  //! fractions
488 };
489 
490 }
491 
492 #endif
virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) const
Get the array of temperature derivatives of the log activity coefficients.
RedlichKisterVPSSTP is a derived class of GibbsExcessVPSSTP that employs the Redlich-Kister approxima...
size_t numBinaryInteractions_
number of binary interaction expressions
virtual void getPartialMolarCp(doublereal *cpbar) const
Returns an array of partial molar entropies 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...
int formRedlichKister_
form of the RedlichKister interaction expression.
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
virtual ThermoPhase * duplMyselfAsThermoPhase() const
Duplication routine for objects which inherit from ThermoPhase.
int formTempModel_
form of the temperature dependence of the Redlich-Kister interaction expression.
Array2D dlnActCoeff_dX_
Two dimensional array of derivatives of activity coefficients wrt mole fractions. ...
virtual std::string type() const
String indicating the thermodynamic model implemented.
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:31
std::vector< size_t > m_pSpecies_A_ij
vector of species indices representing species A in the interaction
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual void getLnActivityCoefficients(doublereal *lnac) const
Get the array of non-dimensional molar-based ln activity coefficients at the current solution tempera...
virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const
Get the array of log species mole number derivatives of the log activity coefficients.
void s_update_dlnActCoeff_dX_() const
Internal routine that calculates the derivative of the activity coefficients wrt the mole fractions...
virtual void getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const
Get the array of temperature second derivatives of the log activity coefficients. ...
void s_update_dlnActCoeff_dT() const
Update the derivative of the log of the activity coefficients wrt T.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
std::vector< size_t > m_pSpecies_B_ij
vector of species indices representing species B in the interaction
Header for intermediate ThermoPhase object for phases which employ Gibbs excess free energy based for...
void readXMLBinarySpecies(XML_Node &xmlBinarySpecies)
Process an XML node called "binaryNeutralSpeciesParameters".
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, etc) along a line in parameter space or along a line in physical space.
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 doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
void initLengths()
Initialize lengths of local variables after all species have been identified.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes 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.
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 Vint(double &VintOut, double &voltsOut)
Utility routine that calculates a literature expression.
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.
Namespace for the Cantera kernel.
Definition: application.cpp:29
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
void resizeNumInteractions(const size_t num)
Resize internal arrays within the object that depend upon the number of binary Redlich-Kister interac...
std::vector< size_t > m_N_ij
Vector of the length of the polynomial for the interaction.
std::vector< vector_fp > m_SE_m_ij
Entropy term for the binary mole fraction interaction of the excess Gibbs free energy expression...
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.