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