diff --git a/Cantera/src/ConstCpPoly.cpp b/Cantera/src/ConstCpPoly.cpp index d0d998385..38e1bcd3e 100644 --- a/Cantera/src/ConstCpPoly.cpp +++ b/Cantera/src/ConstCpPoly.cpp @@ -1,5 +1,5 @@ /** - * @file ConstCpPoly.h + * @file ConstCpPoly.cpp * * $Author$ * $Revision$ @@ -14,119 +14,125 @@ namespace Cantera { - ConstCpPoly::ConstCpPoly() - : m_t0(0.0), - m_cp0_R(0.0), - m_h0_R(0.0), - m_s0_R(0.0), - m_logt0(0.0), - m_lowT(0.0), - m_highT(0.0), - m_Pref(0.0), - m_index(0) { - } + ConstCpPoly::ConstCpPoly() + : m_t0(0.0), + m_cp0_R(0.0), + m_h0_R(0.0), + m_s0_R(0.0), + m_logt0(0.0), + m_lowT(0.0), + m_highT(0.0), + m_Pref(0.0), + m_index(0) { + } - ConstCpPoly::ConstCpPoly(int n, doublereal tlow, doublereal thigh, - doublereal pref, - const doublereal* coeffs) : - m_lowT (tlow), - m_highT (thigh), - m_Pref (pref), - m_index (n) { - m_t0 = coeffs[0]; - m_h0_R = coeffs[1] / GasConstant; - m_s0_R = coeffs[2] / GasConstant; - m_cp0_R = coeffs[3] / GasConstant; - m_logt0 = log(m_t0); - } + ConstCpPoly::ConstCpPoly(int n, doublereal tlow, doublereal thigh, + doublereal pref, + const doublereal* coeffs) : + m_lowT (tlow), + m_highT (thigh), + m_Pref (pref), + m_index (n) { + m_t0 = coeffs[0]; + m_h0_R = coeffs[1] / GasConstant; + m_s0_R = coeffs[2] / GasConstant; + m_cp0_R = coeffs[3] / GasConstant; + m_logt0 = log(m_t0); + } - ConstCpPoly::ConstCpPoly(const ConstCpPoly& b) : - m_t0 (b.m_t0), - m_cp0_R (b.m_cp0_R), - m_h0_R (b.m_h0_R), - m_s0_R (b.m_s0_R), - m_logt0 (b.m_logt0), - m_lowT (b.m_lowT), - m_highT (b.m_highT), - m_Pref (b.m_Pref), - m_index (b.m_index) - { - } + ConstCpPoly::ConstCpPoly(const ConstCpPoly& b) : + m_t0 (b.m_t0), + m_cp0_R (b.m_cp0_R), + m_h0_R (b.m_h0_R), + m_s0_R (b.m_s0_R), + m_logt0 (b.m_logt0), + m_lowT (b.m_lowT), + m_highT (b.m_highT), + m_Pref (b.m_Pref), + m_index (b.m_index) + { + } - ConstCpPoly& ConstCpPoly::operator=(const ConstCpPoly& b) { - if (&b != this) { - m_t0 = b.m_t0; - m_cp0_R = b.m_cp0_R; - m_h0_R = b.m_h0_R; - m_s0_R = b.m_s0_R; - m_logt0 = b.m_logt0; - m_lowT = b.m_lowT; - m_highT = b.m_highT; - m_Pref = b.m_Pref; - m_index = b.m_index; - } - return *this; + ConstCpPoly& ConstCpPoly::operator=(const ConstCpPoly& b) { + if (&b != this) { + m_t0 = b.m_t0; + m_cp0_R = b.m_cp0_R; + m_h0_R = b.m_h0_R; + m_s0_R = b.m_s0_R; + m_logt0 = b.m_logt0; + m_lowT = b.m_lowT; + m_highT = b.m_highT; + m_Pref = b.m_Pref; + m_index = b.m_index; } + return *this; + } - ConstCpPoly::~ConstCpPoly(){} + ConstCpPoly::~ConstCpPoly(){} - SpeciesThermoInterpType * - ConstCpPoly::duplMyselfAsSpeciesThermoInterpType() const { - ConstCpPoly* newCCP = new ConstCpPoly(*this); - return (SpeciesThermoInterpType*) newCCP; - } + SpeciesThermoInterpType * + ConstCpPoly::duplMyselfAsSpeciesThermoInterpType() const { + ConstCpPoly* newCCP = new ConstCpPoly(*this); + return (SpeciesThermoInterpType*) newCCP; + } - doublereal ConstCpPoly::minTemp() const { - return m_lowT; - } - doublereal ConstCpPoly::maxTemp() const { - return m_highT; - } - doublereal ConstCpPoly::refPressure() const { - return m_Pref; - } + doublereal ConstCpPoly::minTemp() const { + return m_lowT; + } + doublereal ConstCpPoly::maxTemp() const { + return m_highT; + } + doublereal ConstCpPoly::refPressure() const { + return m_Pref; + } - void ConstCpPoly::updateProperties(const doublereal* tt, - doublereal* cp_R, - doublereal* h_RT, - doublereal* s_R) const { - double t = *tt; - doublereal logt = log(t); - doublereal rt = 1.0/t; - cp_R[m_index] = m_cp0_R; - h_RT[m_index] = rt*(m_h0_R + (t - m_t0) * m_cp0_R); - s_R[m_index] = m_s0_R + m_cp0_R * (logt - m_logt0); - } + void ConstCpPoly::updateProperties(const doublereal* tt, + doublereal* cp_R, + doublereal* h_RT, + doublereal* s_R) const { + double t = *tt; + doublereal logt = log(t); + doublereal rt = 1.0/t; + cp_R[m_index] = m_cp0_R; + h_RT[m_index] = rt*(m_h0_R + (t - m_t0) * m_cp0_R); + s_R[m_index] = m_s0_R + m_cp0_R * (logt - m_logt0); + } - void ConstCpPoly::updatePropertiesTemp(const doublereal temp, - doublereal* cp_R, - doublereal* h_RT, - doublereal* s_R) const { - doublereal logt = log(temp); - doublereal rt = 1.0/temp; - cp_R[m_index] = m_cp0_R; - h_RT[m_index] = rt*(m_h0_R + (temp - m_t0) * m_cp0_R); - s_R[m_index] = m_s0_R + m_cp0_R * (logt - m_logt0); - } + void ConstCpPoly::updatePropertiesTemp(const doublereal temp, + doublereal* cp_R, + doublereal* h_RT, + doublereal* s_R) const { + doublereal logt = log(temp); + doublereal rt = 1.0/temp; + cp_R[m_index] = m_cp0_R; + h_RT[m_index] = rt*(m_h0_R + (temp - m_t0) * m_cp0_R); + s_R[m_index] = m_s0_R + m_cp0_R * (logt - m_logt0); + } + + void ConstCpPoly::reportParameters(int &n, int &type, + doublereal &tlow, doublereal &thigh, + doublereal &pref, + doublereal* const coeffs) const { + n = m_index; + type = CONSTANT_CP; + tlow = m_lowT; + thigh = m_highT; + pref = m_Pref; + coeffs[0] = m_t0; + coeffs[1] = m_h0_R * GasConstant; + coeffs[2] = m_s0_R * GasConstant; + coeffs[3] = m_cp0_R * GasConstant; + } + + void ConstCpPoly::modifyParameters(doublereal* coeffs) { + m_t0 = coeffs[0]; + m_h0_R = coeffs[1] / GasConstant; + m_s0_R = coeffs[2] / GasConstant; + m_cp0_R = coeffs[3] / GasConstant; + m_logt0 = log(m_t0); + } - void ConstCpPoly::reportParameters(int &n, int &type, - doublereal &tlow, doublereal &thigh, - doublereal &pref, - doublereal* const coeffs) const { - n = m_index; - type = CONSTANT_CP; - tlow = m_lowT; - thigh = m_highT; - pref = m_Pref; - coeffs[0] = m_t0; - coeffs[1] = m_h0_R * GasConstant; - coeffs[2] = m_s0_R * GasConstant; - coeffs[3] = m_cp0_R * GasConstant; - } } - - - diff --git a/Cantera/src/ConstCpPoly.h b/Cantera/src/ConstCpPoly.h index 233901465..7d05ebbc8 100644 --- a/Cantera/src/ConstCpPoly.h +++ b/Cantera/src/ConstCpPoly.h @@ -1,6 +1,9 @@ /** * @file ConstCpPoly.h * + * Declarations for the SpeciesThermoInterpType object that + * employs a constant heat capacity assumption. + * * $Author$ * $Revision$ * $Date$ @@ -16,50 +19,177 @@ namespace Cantera { + /** + * A constant-heat capacity species thermodynamic property manager class. + * This makes the + * assumption that the heat capacity is a constant. Then, the following + * relations are used to complete the specification of the thermodynamic + * functions for the species. + * + * \f[ + * \frac{c_p(T)}{R} = Cp0\_R + * \f] + * \f[ + * \frac{h^0(T)}{RT} = \frac{1}{T} * (h0\_R + (T - T_0) * Cp0\_R) + * \f] + * \f[ + * \frac{s^0(T)}{R} = (s0\_R + (log(T) - log(T_0)) * Cp0\_R) + * \f] + * + * This parameterization takes 4 input values. These are: + * - c[0] = \f$ T_0 \f$(Kelvin) + * - c[1] = \f$ H_k^o(T_0, p_{ref}) \f$ (J/kmol) + * - c[2] = \f$ S_k^o(T_0, p_{ref}) \f$ (J/kmol K) + * - c[3] = \f$ {Cp}_k^o(T_0, p_{ref}) \f$ (J(kmol K) + * + * The multispecies SimpleThermo class makes the same assumptions as + * this class does. + * + * @see SimpleThermo + * + * @ingroup spthermo + */ + class ConstCpPoly: public SpeciesThermoInterpType { - class ConstCpPoly: public SpeciesThermoInterpType { + public: - public: + //! empty constructor + ConstCpPoly(); - ConstCpPoly(); - ConstCpPoly(int n, doublereal tlow, doublereal thigh, - doublereal pref, - const doublereal* coeffs); - ConstCpPoly(const ConstCpPoly&); - ConstCpPoly& operator=(const ConstCpPoly&); - virtual ~ConstCpPoly(); - virtual SpeciesThermoInterpType * - duplMyselfAsSpeciesThermoInterpType() const; + //! Constructor used in templated instantiations + /*! + * @param n Species index + * @param tlow Minimum temperature + * @param thigh Maximum temperature + * @param pref reference pressure (Pa). + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state for species n. + * There are 4 coefficients for the %ConstCpPoly parameterization. + * - c[0] = \f$ T_0 \f$(Kelvin) + * - c[1] = \f$ H_k^o(T_0, p_{ref}) \f$ (J/kmol) + * - c[2] = \f$ S_k^o(T_0, p_{ref}) \f$ (J/kmol K) + * - c[3] = \f$ {Cp}_k^o(T_0, p_{ref}) \f$ (J(kmol K) + * + */ + ConstCpPoly(int n, doublereal tlow, doublereal thigh, + doublereal pref, + const doublereal* coeffs); - doublereal minTemp() const; - doublereal maxTemp() const; - doublereal refPressure() const; - virtual int reportType() const { return CONSTANT_CP; } - - void updateProperties(const doublereal* tt, - doublereal* cp_R, doublereal* h_RT, + //! copy constructor + ConstCpPoly(const ConstCpPoly&); + + //! Assignment operator + ConstCpPoly& operator=(const ConstCpPoly&); + + //! Destructor + virtual ~ConstCpPoly(); + + //! Duplicator + virtual SpeciesThermoInterpType * + duplMyselfAsSpeciesThermoInterpType() const; + //! Returns the minimum temperature that the thermo + //! parameterization is valid + doublereal minTemp() const; + + //! Returns the maximum temperature that the thermo + //! parameterization is valid + doublereal maxTemp() const; + + //! Returns the reference pressure (Pa) + doublereal refPressure() const; + + //! Returns an integer representing the type of parameterization + virtual int reportType() const { return CONSTANT_CP; } + + //! Update the properties for this species, given a temperature polynomial + /*! + * This method is called with a pointer to an array containing the functions of + * temperature needed by this parameterization, and three pointers to arrays where the + * computed property values should be written. This method updates only one value in + * each array. + * + * Form and Length of the temperature polynomial: + * - m_t[0] = tt; + * + * @param tt Vector of temperature polynomials + * @param cp_R Vector of Dimensionless heat capacities. + * (length m_kk). + * @param h_RT Vector of Dimensionless enthalpies. + * (length m_kk). + * @param s_R Vector of Dimensionless entropies. + * (length m_kk). + */ + void updateProperties(const doublereal* tt, + doublereal* cp_R, doublereal* h_RT, + doublereal* s_R) const; + + //! Compute the reference-state property of one species + /*! + * Given temperature T in K, this method updates the values of + * the non-dimensional heat capacity at constant pressure, + * enthalpy, and entropy, at the reference pressure, Pref + * of one of the species. The species index is used + * to reference into the cp_R, h_RT, and s_R arrays. + * + * @param temp Temperature (Kelvin) + * @param cp_R Vector of Dimensionless heat capacities. + * (length m_kk). + * @param h_RT Vector of Dimensionless enthalpies. + * (length m_kk). + * @param s_R Vector of Dimensionless entropies. + * (length m_kk). + */ + void updatePropertiesTemp(const doublereal temp, + doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const; + //!This utility function reports back the type of + //! parameterization and all of the parameters for the + //! species, index. + /*! + * All parameters are output variables + * + * @param n Species index + * @param type Integer type of the standard type + * @param tlow output - Minimum temperature + * @param thigh output - Maximum temperature + * @param pref output - reference pressure (Pa). + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state. + */ + void reportParameters(int &n, int &type, + doublereal &tlow, doublereal &thigh, + doublereal &pref, + doublereal* const coeffs) const; + //! Modify parameters for the standard state + /*! + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state. + */ + virtual void modifyParameters(doublereal* coeffs); - void updatePropertiesTemp(const doublereal temp, - doublereal* cp_R, doublereal* h_RT, - doublereal* s_R) const; - - void reportParameters(int &n, int &type, - doublereal &tlow, doublereal &thigh, - doublereal &pref, - doublereal* const coeffs) const; - protected: - doublereal m_t0; - doublereal m_cp0_R; - doublereal m_h0_R; - doublereal m_s0_R; - doublereal m_logt0; - doublereal m_lowT, m_highT, m_Pref; - int m_index; + protected: + //! Base temperature + doublereal m_t0; + //! Dimensionless value of the heat capacity + doublereal m_cp0_R; + //! dimensionless value of the enthaply at t0 + doublereal m_h0_R; + //! Dimensionless value of the entropy at t0 + doublereal m_s0_R; + //! log of the t0 value + doublereal m_logt0; + //! Minimum temperature for which the parameterization is valid (Kelvin) + doublereal m_lowT; + //! Maximum temperature for which the parameterization is valid (Kelvin) + doublereal m_highT; + //! Reference pressure (Pa) + doublereal m_Pref; + //! Species Index + int m_index; - private: + private: - }; + }; } diff --git a/Cantera/src/GeneralSpeciesThermo.cpp b/Cantera/src/GeneralSpeciesThermo.cpp index a5fa0e3b5..83490f3ee 100644 --- a/Cantera/src/GeneralSpeciesThermo.cpp +++ b/Cantera/src/GeneralSpeciesThermo.cpp @@ -69,7 +69,7 @@ namespace Cantera { } - /** + /* * Install parameterization for a species. * @param index Species index * @param type parameterization type diff --git a/Cantera/src/Mu0Poly.cpp b/Cantera/src/Mu0Poly.cpp index ea33c0648..a779668c5 100644 --- a/Cantera/src/Mu0Poly.cpp +++ b/Cantera/src/Mu0Poly.cpp @@ -1,5 +1,7 @@ /** - * @file Mu0Poly.h + * @file Mu0Poly.cpp + * + * Definitions for the Mu0Poly class. * * $Author$ * $Revision$ @@ -19,389 +21,388 @@ using namespace ctml; namespace Cantera { - /** - * The Mu0Poly class implements a linear interpolation - * of the standard state chemical potential of one - * species at a single reference pressure. - * The chemical potential is input as a series of (T, mu0) - * values. The first temperature is assumed to be equal - * to 298.15 K; however, this may be relaxed in the future. - * This information, and an assumption of a constant - * heat capacity within each interval is enough to - * calculate all thermodynamic functions. - * - * The basic equation for going from point 1 to point 2 - * are as follows for T, T1 <= T <= T2 - * - * mu1 = H1 - T1 * S1 - * - * mu2 - mu1 = Cp1(T2 - T1) - Cp1(ln(T2/T1)) - S1(T2 - T1) - * - * S2 = S1 + Cp1(ln(T2/T1)) - * - * H2 = H1 + Cp1(T2 - T1) - * - * In the future, a better assumption about the heat - * capacity may be employed, so that it can be continuous. - * - * Notes about temperature interpolation for T < T1 and T > Tn - * These are achieved by assuming a constant heat capacity - * equal to the value in the closest temperature interval. - * No error is thrown. - */ - Mu0Poly::Mu0Poly() : m_numIntervals(0), - m_H298(0.0), - m_lowT(0.0), - m_highT(0.0), - m_Pref(0.0), - m_index(0) { - } - - /** - * Mu0Poly(): - * - * In the constructor, we calculate and store the - * piecewise linear approximation to the thermodynamic - * functions. - * - * coeffs[0] = number of points (integer) - * 1 = H298(J/kmol) - * 2 = T1 (Kelvin) - * 3 = mu1 (J/kmol) - * 4 = T2 (Kelvin) - * 5 = mu2 (J/kmol) - * 6 = T3 (Kelvin) - * 7 = mu3 (J/kmol) - * ........ - */ - Mu0Poly::Mu0Poly(int n, doublereal tlow, doublereal thigh, - doublereal pref, - const doublereal* coeffs) : - m_numIntervals(0), - m_H298(0.0), - m_lowT (tlow), - m_highT (thigh), - m_Pref (pref), - m_index (n) { - - int i, iindex; - double T1, T2; - int nPoints = (int) coeffs[0]; - if (nPoints < 2) { - throw CanteraError("Mu0Poly", - "nPoints must be >= 2"); - } - m_numIntervals = nPoints - 1; - m_H298 = coeffs[1] / GasConstant; - int iT298 = 0; - /* - * Resize according to the number of points - */ - m_t0_int.resize(nPoints); - m_h0_R_int.resize(nPoints); - m_s0_R_int.resize(nPoints); - m_cp0_R_int.resize(nPoints); - m_mu0_R_int.resize(nPoints); - /* - * Calculate the T298 interval and make sure that - * the temperatures are strictly monotonic. - * Also distribute the data into the internal arrays. - */ - bool ifound = false; - for (i = 0, iindex = 2; i < nPoints; i++) { - T1 = coeffs[iindex]; - m_t0_int[i] = T1; - m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant; - if (T1 == 298.15) { - iT298 = i; - ifound = true; - } - if (i < nPoints - 1) { - T2 = coeffs[iindex+2]; - if (T2 <= T1) { - throw CanteraError("Mu0Poly", - "Temperatures are not monotonic increasing"); - } - } - iindex += 2; - } - if (!ifound) { - throw CanteraError("Mu0Poly", - "One temperature has to be 298.15"); - } - - /* - * Starting from the interval with T298, we go up - */ - doublereal mu2, s1, s2, h1, h2, cpi, deltaMu, deltaT; - T1 = m_t0_int[iT298]; - doublereal mu1 = m_mu0_R_int[iT298]; - m_h0_R_int[iT298] = m_H298; - m_s0_R_int[iT298] = - (mu1 - m_h0_R_int[iT298]) / T1; - for (i = iT298; i < m_numIntervals; i++) { - T1 = m_t0_int[i]; - s1 = m_s0_R_int[i]; - h1 = m_h0_R_int[i]; - mu1 = m_mu0_R_int[i]; - T2 = m_t0_int[i+1]; - mu2 = m_mu0_R_int[i+1]; - deltaMu = mu2 - mu1; - deltaT = T2 - T1; - cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1)); - h2 = h1 + cpi * deltaT; - s2 = s1 + cpi * log(T2/T1); - m_cp0_R_int[i] = cpi; - m_h0_R_int[i+1] = h2; - m_s0_R_int[i+1] = s2; - m_cp0_R_int[i+1] = cpi; - } - - /* - * Starting from the interval with T298, we go down - */ - if (iT298 > 0) { - T2 = m_t0_int[iT298]; - mu2 = m_mu0_R_int[iT298]; - m_h0_R_int[iT298] = m_H298; - m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2; - for (i = iT298 - 1; i >= 0; i--) { - T1 = m_t0_int[i]; - mu1 = m_mu0_R_int[i]; - T2 = m_t0_int[i+1]; - mu2 = m_mu0_R_int[i+1]; - s2 = m_s0_R_int[i+1]; - h2 = m_h0_R_int[i+1]; - deltaMu = mu2 - mu1; - deltaT = T2 - T1; - cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1)); - h1 = h2 - cpi * deltaT; - s1 = s2 - cpi * log(T2/T1); - m_cp0_R_int[i] = cpi; - m_h0_R_int[i] = h1; - m_s0_R_int[i] = s1; - if (i == (m_numIntervals-1)) { - m_cp0_R_int[i+1] = cpi; - } - } - } -#ifdef DEBUG_HKM_NOT - printf(" Temp mu0(J/kmol) cp0(J/kmol/K) " - " h0(J/kmol) s0(J/kmol/K) \n"); - for (i = 0; i < nPoints; i++) { - printf("%12.3g %12.5g %12.5g %12.5g %12.5g\n", - m_t0_int[i], m_mu0_R_int[i] * GasConstant, - m_cp0_R_int[i]* GasConstant, - m_h0_R_int[i]* GasConstant, - m_s0_R_int[i]* GasConstant); - fflush(stdout); - } -#endif - } - - - Mu0Poly::Mu0Poly(const Mu0Poly &b) - : m_numIntervals (b.m_numIntervals), - m_H298 (b.m_H298), - m_t0_int (b.m_t0_int), - m_mu0_R_int (b.m_mu0_R_int), - m_h0_R_int (b.m_h0_R_int), - m_s0_R_int (b.m_s0_R_int), - m_cp0_R_int (b.m_cp0_R_int), - m_lowT (b.m_lowT), - m_highT (b.m_highT), - m_Pref (b.m_Pref), - m_index (b.m_index) { - } - - Mu0Poly& Mu0Poly::operator=(const Mu0Poly& b) { - if (&b != this) { - m_numIntervals = b.m_numIntervals; - m_H298 = b.m_H298; - m_t0_int = b.m_t0_int; - m_mu0_R_int = b.m_mu0_R_int; - m_h0_R_int = b.m_h0_R_int; - m_s0_R_int = b.m_s0_R_int; - m_cp0_R_int = b.m_cp0_R_int; - m_lowT = b.m_lowT; - m_highT = b.m_highT; - m_Pref = b.m_Pref; - m_index = b.m_index; - } - return *this; + + Mu0Poly::Mu0Poly() : m_numIntervals(0), + m_H298(0.0), + m_lowT(0.0), + m_highT(0.0), + m_Pref(0.0), + m_index(0) { + } + + /* + * Mu0Poly(): + * + * In the constructor, we calculate and store the + * piecewise linear approximation to the thermodynamic + * functions. + * + * coeffs[0] = number of points (integer) + * 1 = H298(J/kmol) + * 2 = T1 (Kelvin) + * 3 = mu1 (J/kmol) + * 4 = T2 (Kelvin) + * 5 = mu2 (J/kmol) + * 6 = T3 (Kelvin) + * 7 = mu3 (J/kmol) + * ........ + */ + Mu0Poly::Mu0Poly(int n, doublereal tlow, doublereal thigh, + doublereal pref, + const doublereal* coeffs) : + m_numIntervals(0), + m_H298(0.0), + m_lowT (tlow), + m_highT (thigh), + m_Pref (pref), + m_index (n) { + + processCoeffs(coeffs); + } + + + Mu0Poly::Mu0Poly(const Mu0Poly &b) + : m_numIntervals (b.m_numIntervals), + m_H298 (b.m_H298), + m_t0_int (b.m_t0_int), + m_mu0_R_int (b.m_mu0_R_int), + m_h0_R_int (b.m_h0_R_int), + m_s0_R_int (b.m_s0_R_int), + m_cp0_R_int (b.m_cp0_R_int), + m_lowT (b.m_lowT), + m_highT (b.m_highT), + m_Pref (b.m_Pref), + m_index (b.m_index) { + } + + Mu0Poly& Mu0Poly::operator=(const Mu0Poly& b) { + if (&b != this) { + m_numIntervals = b.m_numIntervals; + m_H298 = b.m_H298; + m_t0_int = b.m_t0_int; + m_mu0_R_int = b.m_mu0_R_int; + m_h0_R_int = b.m_h0_R_int; + m_s0_R_int = b.m_s0_R_int; + m_cp0_R_int = b.m_cp0_R_int; + m_lowT = b.m_lowT; + m_highT = b.m_highT; + m_Pref = b.m_Pref; + m_index = b.m_index; } + return *this; + } - /** - * Destructor: - */ - Mu0Poly::~Mu0Poly(){ - } + /** + * Destructor: + */ + Mu0Poly::~Mu0Poly(){ + } - SpeciesThermoInterpType * - Mu0Poly::duplMyselfAsSpeciesThermoInterpType() const { - Mu0Poly* mp = new Mu0Poly(*this); - return (SpeciesThermoInterpType *) mp; - } + SpeciesThermoInterpType * + Mu0Poly::duplMyselfAsSpeciesThermoInterpType() const { + Mu0Poly* mp = new Mu0Poly(*this); + return (SpeciesThermoInterpType *) mp; + } - doublereal Mu0Poly::minTemp() const { return m_lowT;} - doublereal Mu0Poly::maxTemp() const { return m_highT;} - doublereal Mu0Poly::refPressure() const { return m_Pref; } + doublereal Mu0Poly::minTemp() const { return m_lowT;} + doublereal Mu0Poly::maxTemp() const { return m_highT;} + doublereal Mu0Poly::refPressure() const { return m_Pref; } - /** - * updateProperties is the main workhorse program. - * Given a temperature (*tt), it calculates the thermodynamic - * functions H/RT, S_R, and cp_R, and returns the answer. - * - * Note, it returns an answer by inserting the values into the - * index position, m_index in vectors of H/RT, S_R, and cp_R. - * - * - * Input - * ------- - * *tt = Temperature (Kelvin) - * - */ - void Mu0Poly:: - updateProperties(const doublereal* tt, doublereal* cp_R, - doublereal* h_RT, doublereal* s_R) const { - int j = m_numIntervals; - double T = *tt; - for (int i = 0; i < m_numIntervals; i++) { - double T2 = m_t0_int[i+1]; - if (T <=T2) { - j = i; - break; - } - } - double T1 = m_t0_int[j]; - double cp_Rj = m_cp0_R_int[j]; - - doublereal rt = 1.0/T; - cp_R[m_index] = cp_Rj; - h_RT[m_index] = rt*(m_h0_R_int[j] + (T - T1) * cp_Rj); - s_R[m_index] = m_s0_R_int[j] + cp_Rj * (log(T/T1)); + /** + * updateProperties is the main workhorse program. + * Given a temperature (*tt), it calculates the thermodynamic + * functions H/RT, S_R, and cp_R, and returns the answer. + * + * Note, it returns an answer by inserting the values into the + * index position, m_index in vectors of H/RT, S_R, and cp_R. + * + * + * Input + * ------- + * *tt = Temperature (Kelvin) + * + */ + void Mu0Poly:: + updateProperties(const doublereal* tt, doublereal* cp_R, + doublereal* h_RT, doublereal* s_R) const { + int j = m_numIntervals; + double T = *tt; + for (int i = 0; i < m_numIntervals; i++) { + double T2 = m_t0_int[i+1]; + if (T <=T2) { + j = i; + break; + } } + double T1 = m_t0_int[j]; + double cp_Rj = m_cp0_R_int[j]; - void Mu0Poly:: - updatePropertiesTemp(const doublereal T, - doublereal* cp_R, - doublereal* h_RT, - doublereal* s_R) const { - updateProperties(&T, cp_R, h_RT, s_R); + doublereal rt = 1.0/T; + cp_R[m_index] = cp_Rj; + h_RT[m_index] = rt*(m_h0_R_int[j] + (T - T1) * cp_Rj); + s_R[m_index] = m_s0_R_int[j] + cp_Rj * (log(T/T1)); + } + + void Mu0Poly:: + updatePropertiesTemp(const doublereal T, + doublereal* cp_R, + doublereal* h_RT, + doublereal* s_R) const { + updateProperties(&T, cp_R, h_RT, s_R); + } + + /* + * report all of the parameters that make up this + * interpolation. + * + * + */ + void Mu0Poly::reportParameters(int &n, int &type, + doublereal &tlow, doublereal &thigh, + doublereal &pref, + doublereal* const coeffs) const { + n = m_index; + type = MU0_INTERP; + tlow = m_lowT; + thigh = m_highT; + pref = m_Pref; + coeffs[0] = m_numIntervals+1; + coeffs[1] = m_H298 * GasConstant; + int j = 2; + for (int i = 0; i < m_numIntervals+1; i++) { + coeffs[j] = m_t0_int[i]; + coeffs[j+1] = m_mu0_R_int[i] * GasConstant; + j += 2; } + } - /** - * report all of the parameters that make up this - * interpolation. - */ - void Mu0Poly::reportParameters(int &n, int &type, - doublereal &tlow, doublereal &thigh, - doublereal &pref, - doublereal* const coeffs) const { - n = m_index; - type = MU0_INTERP; - tlow = m_lowT; - thigh = m_highT; - pref = m_Pref; - coeffs[0] = m_numIntervals+1; - coeffs[1] = m_H298 * GasConstant; - int j = 2; - for (int i = 0; i < m_numIntervals+1; i++) { - coeffs[j] = m_t0_int[i]; - coeffs[j+1] = m_mu0_R_int[i] * GasConstant; - j += 2; - } - } + void Mu0Poly::modifyParameters(doublereal* coeffs) { + processCoeffs(coeffs); + } + /* + * Install a Mu0 polynomial thermodynamic reference state property + * parameterization for species k into a SpeciesThermo instance, + * getting the information from an XML database. + */ + void installMu0ThermoFromXML(std::string speciesName, + SpeciesThermo& sp, int k, + const XML_Node* Mu0Node_ptr) { - /** - * Install a Mu0 polynomial thermodynamic reference state property - * parameterization for species k into a SpeciesThermo instance, - * getting the information from an XML database. - */ - void installMu0ThermoFromXML(std::string speciesName, - SpeciesThermo& sp, int k, - const XML_Node* Mu0Node_ptr) { - - doublereal tmin, tmax; - bool dimensionlessMu0Values = false; - const XML_Node& Mu0Node = *Mu0Node_ptr; + doublereal tmin, tmax; + bool dimensionlessMu0Values = false; + const XML_Node& Mu0Node = *Mu0Node_ptr; - tmin = fpValue(Mu0Node["Tmin"]); - tmax = fpValue(Mu0Node["Tmax"]); - doublereal pref = fpValue(Mu0Node["Pref"]); + tmin = fpValue(Mu0Node["Tmin"]); + tmax = fpValue(Mu0Node["Tmax"]); + doublereal pref = fpValue(Mu0Node["Pref"]); - doublereal h298 = 0.0; - if (Mu0Node.hasChild("H298")) { - h298 = getFloat(Mu0Node, "H298", "actEnergy"); - } - - int numPoints = 1; - if (Mu0Node.hasChild("numPoints")) { - numPoints = getInteger(Mu0Node, "numPoints"); - } - - vector_fp cValues(numPoints); - const XML_Node *valNode_ptr = - getByTitle(const_cast(Mu0Node), "Mu0Values"); - if (!valNode_ptr) { - throw CanteraError("installMu0ThermoFromXML", - "missing required while processing " - + speciesName); - } - getFloatArray(*valNode_ptr, cValues, true, "actEnergy"); - /* - * Check to see whether the Mu0's were input in a dimensionless - * form. If they were, then the assumed temperature needs to be - * adjusted from the assumed T = 273.15 - */ - string uuu = (*valNode_ptr)["units"]; - if (uuu == "Dimensionless") { - dimensionlessMu0Values = true; - } - int ns = cValues.size(); - if (ns != numPoints) { - throw CanteraError("installMu0ThermoFromXML", - "numPoints inconsistent while processing " - + speciesName); - } - - vector_fp cTemperatures(numPoints); - const XML_Node *tempNode_ptr = - getByTitle(const_cast(Mu0Node), "Mu0Temperatures"); - if (!tempNode_ptr) { - throw CanteraError("installMu0ThermoFromXML", - "missing required while processing + " - + speciesName); - } - getFloatArray(*tempNode_ptr, cTemperatures, false); - ns = cTemperatures.size(); - if (ns != numPoints) { - throw CanteraError("installMu0ThermoFromXML", - "numPoints inconsistent while processing " - + speciesName); - } - - /* - * Fix up dimensionless Mu0 values if input - */ - if (dimensionlessMu0Values) { - for (int i = 0; i < numPoints; i++) { - cValues[i] *= cTemperatures[i] / 273.15; - } - } - - - vector_fp c(2 + 2 * numPoints); - - c[0] = numPoints; - c[1] = h298; - for (int i = 0; i < numPoints; i++) { - c[2+i*2] = cTemperatures[i]; - c[2+i*2+1] = cValues[i]; - } - - sp.install(speciesName, k, MU0_INTERP, &c[0], tmin, tmax, pref); + doublereal h298 = 0.0; + if (Mu0Node.hasChild("H298")) { + h298 = getFloat(Mu0Node, "H298", "actEnergy"); } + + int numPoints = 1; + if (Mu0Node.hasChild("numPoints")) { + numPoints = getInteger(Mu0Node, "numPoints"); + } + + vector_fp cValues(numPoints); + const XML_Node *valNode_ptr = + getByTitle(const_cast(Mu0Node), "Mu0Values"); + if (!valNode_ptr) { + throw CanteraError("installMu0ThermoFromXML", + "missing required while processing " + + speciesName); + } + getFloatArray(*valNode_ptr, cValues, true, "actEnergy"); + /* + * Check to see whether the Mu0's were input in a dimensionless + * form. If they were, then the assumed temperature needs to be + * adjusted from the assumed T = 273.15 + */ + string uuu = (*valNode_ptr)["units"]; + if (uuu == "Dimensionless") { + dimensionlessMu0Values = true; + } + int ns = cValues.size(); + if (ns != numPoints) { + throw CanteraError("installMu0ThermoFromXML", + "numPoints inconsistent while processing " + + speciesName); + } + + vector_fp cTemperatures(numPoints); + const XML_Node *tempNode_ptr = + getByTitle(const_cast(Mu0Node), "Mu0Temperatures"); + if (!tempNode_ptr) { + throw CanteraError("installMu0ThermoFromXML", + "missing required while processing + " + + speciesName); + } + getFloatArray(*tempNode_ptr, cTemperatures, false); + ns = cTemperatures.size(); + if (ns != numPoints) { + throw CanteraError("installMu0ThermoFromXML", + "numPoints inconsistent while processing " + + speciesName); + } + + /* + * Fix up dimensionless Mu0 values if input + */ + if (dimensionlessMu0Values) { + for (int i = 0; i < numPoints; i++) { + cValues[i] *= cTemperatures[i] / 273.15; + } + } + + + vector_fp c(2 + 2 * numPoints); + + c[0] = numPoints; + c[1] = h298; + for (int i = 0; i < numPoints; i++) { + c[2+i*2] = cTemperatures[i]; + c[2+i*2+1] = cValues[i]; + } + + sp.install(speciesName, k, MU0_INTERP, &c[0], tmin, tmax, pref); + } + + /* + * Mu0Poly(): + * + * In the constructor, we calculate and store the + * piecewise linear approximation to the thermodynamic + * functions. + * + * coeffs[0] = number of points (integer) + * 1 = H298(J/kmol) + * 2 = T1 (Kelvin) + * 3 = mu1 (J/kmol) + * 4 = T2 (Kelvin) + * 5 = mu2 (J/kmol) + * 6 = T3 (Kelvin) + * 7 = mu3 (J/kmol) + * ........ + */ + void Mu0Poly::processCoeffs(const doublereal* coeffs) { + + int i, iindex; + double T1, T2; + int nPoints = (int) coeffs[0]; + if (nPoints < 2) { + throw CanteraError("Mu0Poly", + "nPoints must be >= 2"); + } + m_numIntervals = nPoints - 1; + m_H298 = coeffs[1] / GasConstant; + int iT298 = 0; + /* + * Resize according to the number of points + */ + m_t0_int.resize(nPoints); + m_h0_R_int.resize(nPoints); + m_s0_R_int.resize(nPoints); + m_cp0_R_int.resize(nPoints); + m_mu0_R_int.resize(nPoints); + /* + * Calculate the T298 interval and make sure that + * the temperatures are strictly monotonic. + * Also distribute the data into the internal arrays. + */ + bool ifound = false; + for (i = 0, iindex = 2; i < nPoints; i++) { + T1 = coeffs[iindex]; + m_t0_int[i] = T1; + m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant; + if (T1 == 298.15) { + iT298 = i; + ifound = true; + } + if (i < nPoints - 1) { + T2 = coeffs[iindex+2]; + if (T2 <= T1) { + throw CanteraError("Mu0Poly", + "Temperatures are not monotonic increasing"); + } + } + iindex += 2; + } + if (!ifound) { + throw CanteraError("Mu0Poly", + "One temperature has to be 298.15"); + } + + /* + * Starting from the interval with T298, we go up + */ + doublereal mu2, s1, s2, h1, h2, cpi, deltaMu, deltaT; + T1 = m_t0_int[iT298]; + doublereal mu1 = m_mu0_R_int[iT298]; + m_h0_R_int[iT298] = m_H298; + m_s0_R_int[iT298] = - (mu1 - m_h0_R_int[iT298]) / T1; + for (i = iT298; i < m_numIntervals; i++) { + T1 = m_t0_int[i]; + s1 = m_s0_R_int[i]; + h1 = m_h0_R_int[i]; + mu1 = m_mu0_R_int[i]; + T2 = m_t0_int[i+1]; + mu2 = m_mu0_R_int[i+1]; + deltaMu = mu2 - mu1; + deltaT = T2 - T1; + cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1)); + h2 = h1 + cpi * deltaT; + s2 = s1 + cpi * log(T2/T1); + m_cp0_R_int[i] = cpi; + m_h0_R_int[i+1] = h2; + m_s0_R_int[i+1] = s2; + m_cp0_R_int[i+1] = cpi; + } + + /* + * Starting from the interval with T298, we go down + */ + if (iT298 > 0) { + T2 = m_t0_int[iT298]; + mu2 = m_mu0_R_int[iT298]; + m_h0_R_int[iT298] = m_H298; + m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2; + for (i = iT298 - 1; i >= 0; i--) { + T1 = m_t0_int[i]; + mu1 = m_mu0_R_int[i]; + T2 = m_t0_int[i+1]; + mu2 = m_mu0_R_int[i+1]; + s2 = m_s0_R_int[i+1]; + h2 = m_h0_R_int[i+1]; + deltaMu = mu2 - mu1; + deltaT = T2 - T1; + cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1)); + h1 = h2 - cpi * deltaT; + s1 = s2 - cpi * log(T2/T1); + m_cp0_R_int[i] = cpi; + m_h0_R_int[i] = h1; + m_s0_R_int[i] = s1; + if (i == (m_numIntervals-1)) { + m_cp0_R_int[i+1] = cpi; + } + } + } +#ifdef DEBUG_HKM_NOT + printf(" Temp mu0(J/kmol) cp0(J/kmol/K) " + " h0(J/kmol) s0(J/kmol/K) \n"); + for (i = 0; i < nPoints; i++) { + printf("%12.3g %12.5g %12.5g %12.5g %12.5g\n", + m_t0_int[i], m_mu0_R_int[i] * GasConstant, + m_cp0_R_int[i]* GasConstant, + m_h0_R_int[i]* GasConstant, + m_s0_R_int[i]* GasConstant); + fflush(stdout); + } +#endif + } + } diff --git a/Cantera/src/Mu0Poly.h b/Cantera/src/Mu0Poly.h index 756ce9171..e4c87ff23 100644 --- a/Cantera/src/Mu0Poly.h +++ b/Cantera/src/Mu0Poly.h @@ -1,5 +1,7 @@ /** * @file Mu0Poly.h + * + * Declarations for piecewise constant mu0 interpolation. */ /* $Author$ @@ -15,123 +17,289 @@ #include "SpeciesThermoInterpType.h" namespace Cantera { - class SpeciesThermo; - class XML_Node; - /** - * The Mu0Poly class implements a linear interpolation - * of the standard state chemical potential of one - * species at a single reference pressure. - * The chemical potential is input as a series of (T, mu0) - * values. The first temperature is assumed to be equal - * to 298.15 K; however, this may be relaxed in the future. - * This information, and an assumption of a constant - * heat capacity within each interval is enough to - * calculate all thermodynamic functions. - * - * The basic equation for going from point 1 to point 2 - * are as follows for T, T1 <= T <= T2 + class SpeciesThermo; + class XML_Node; + + //! The %Mu0Poly class implements an interpolation of the Gibbs free energy based on a + //! piecewise constant heat capacity approximation. + /*! + * The %Mu0Poly class implements a piecewise constant heat capacity approximation. + * of the standard state chemical potential of one + * species at a single reference pressure. + * The chemical potential is input as a series of (\f$T\f$, \f$ \mu^o(T)\f$) + * values. The first temperature is assumed to be equal + * to 298.15 K; however, this may be relaxed in the future. + * This information, and an assumption of a constant + * heat capacity within each interval is enough to + * calculate all thermodynamic functions. + * + * The piece-wise constant heat capacity is calculated from the change in the chemical potential over each interval. + * Once the heat capacity is known, the other thermodynamic functions may be determined. + * The basic equation for going from temperature point 1 to temperature point 2 + * are as follows for \f$ T \f$, \f$ T_1 <= T <= T_2 \f$ + * + * \f[ + * \mu^o(T_1) = h^o(T_1) - T_1 * s^o(T_1) + * \f] + * \f[ + * \mu^o(T_2) - \mu^o(T_1) = Cp^o(T_1)(T_2 - T_1) - Cp^o(T_1)(T_2)ln(\frac{T_2}{T_1}) - s^o(T_1)(T_2 - T_1) + * \f] + * \f[ + * s^o(T_2) = s^o(T_1) + Cp^o(T_1)ln(\frac{T_2}{T_1}) + * \f] + * \f[ + * h^o(T_2) = h^o(T_1) + Cp^o(T_1)(T_2 - T_1) + * \f] + * + * Within each interval the following relations are used. For \f$ T \f$, \f$ T_1 <= T <= T_2 \f$ + * + * \f[ + * \mu^o(T) = \mu^o(T_1) + Cp^o(T_1)(T - T_1) - Cp^o(T_1)(T_2)ln(\frac{T}{T_1}) - s^o(T_1)(T - T_1) + * \f] + * \f[ + * s^o(T) = s^o(T_1) + Cp^o(T_1)ln(\frac{T}{T_1}) + * \f] + * \f[ + * h^o(T) = h^o(T_1) + Cp^o(T_1)(T - T_1) + * \f] + * + * Notes about temperature interpolation for \f$ T < T_1 \f$ and \f$ T > T_{npoints} \f$. + * These are achieved by assuming a constant heat capacity + * equal to the value in the closest temperature interval. + * No error is thrown. + * + * @note In the future, a better assumption about the heat + * capacity may be employed, so that it can be continuous. + * + * @ingroup spthermo + */ + class Mu0Poly: public SpeciesThermoInterpType { + + public: + + //! Constructor + Mu0Poly(); + + //! Constructor used in templated instantiations + /*! * - * mu1 = H1 - T1 * S1 - * - * mu2 - mu1 = Cp1(T2 - T1) - Cp1(ln(T2/T1)) - S1(T2 - T1) + * In the constructor, we calculate and store the + * piecewise linear approximation to the thermodynamic + * functions. * - * S2 = S1 + Cp1(ln(T2/T1)) - * - * H2 = H1 + Cp1(T2 - T1) - * - * In the future, a better assumption about the heat - * capacity may be employed, so that it can be continuous. - * - * Notes about temperature interpolation for T < T1 and T > Tn - * These are achieved by assuming a constant heat capacity - * equal to the value in the closest temperature interval. - * No error is thrown. + * @param n Species index + * @param tlow Minimum temperature + * @param thigh Maximum temperature + * @param pref reference pressure (Pa). + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state for species n. + * There are \f$ 2+npoints*2 \f$ coefficients, where + * \f$ npoints \f$ are the number of temperature points. + * Their identity is further broken down: + * - coeffs[0] = number of points (integer) + * - coeffs[1] = \f$ h^o(298.15 K) \f$ (J/kmol) + * - coeffs[2] = \f$ T_1 \f$ (Kelvin) + * - coeffs[3] = \f$ \mu^o(T_1) \f$ (J/kmol) + * - coeffs[4] = \f$ T_2 \f$ (Kelvin) + * - coeffs[5] = \f$ \mu^o(T_2) \f$ (J/kmol) + * - coeffs[6] = \f$ T_3 \f$ (Kelvin) + * - coeffs[7] = \f$ \mu^o(T_3) \f$ (J/kmol) + * - ........ + * . */ - class Mu0Poly: public SpeciesThermoInterpType { + Mu0Poly(int n, doublereal tlow, doublereal thigh, + doublereal pref, const doublereal* coeffs); - public: + //! Copy constructor + Mu0Poly(const Mu0Poly &); - Mu0Poly(); + //! Assignment operator + Mu0Poly& operator=(const Mu0Poly&); - Mu0Poly(int n, doublereal tlow, doublereal thigh, - doublereal pref, const doublereal* coeffs); - Mu0Poly(const Mu0Poly &); - Mu0Poly& operator=(const Mu0Poly&); - virtual ~Mu0Poly(); - SpeciesThermoInterpType * - duplMyselfAsSpeciesThermoInterpType() const; - - doublereal minTemp() const; - doublereal maxTemp() const; - doublereal refPressure() const; - virtual int reportType() const { return MU0_INTERP; } + //! Destructor + virtual ~Mu0Poly(); - /** - * Update all of the properties, using the polynomial - * tPoly[] - * - * tPoly[0] = temp (Kelvin) - */ - void updateProperties(const doublereal* tPoly, - doublereal* cp_R, doublereal* h_RT, + //! Duplicator + virtual SpeciesThermoInterpType * + duplMyselfAsSpeciesThermoInterpType() const; + + //! Returns the minimum temperature that the thermo + //! parameterization is valid + virtual doublereal minTemp() const; + + //! Returns the maximum temperature that the thermo + //! parameterization is valid + virtual doublereal maxTemp() const; + + //! Returns the reference pressure (Pa) + virtual doublereal refPressure() const; + + //! Returns an integer representing the type of parameterization + virtual int reportType() const { return MU0_INTERP; } + + + //! Update the properties for this species, given a temperature polynomial + /*! + * This method is called with a pointer to an array containing the functions of + * temperature needed by this parameterization, and three pointers to arrays where the + * computed property values should be written. This method updates only one value in + * each array. + * + * Temperature Polynomial: + * + * tPoly[0] = temp (Kelvin) + * + * @param tPoly vector of temperature polynomials. Length = 1 + * @param cp_R Vector of Dimensionless heat capacities. + * (length m_kk). + * @param h_RT Vector of Dimensionless enthalpies. + * (length m_kk). + * @param s_R Vector of Dimensionless entropies. + * (length m_kk). + */ + virtual void updateProperties(const doublereal* tPoly, + doublereal* cp_R, doublereal* h_RT, + doublereal* s_R) const ; + + //! Compute the reference-state property of one species + /*! + * Given temperature T in K, this method updates the values of + * the non-dimensional heat capacity at constant pressure, + * enthalpy, and entropy, at the reference pressure, Pref + * of one of the species. The species index is used + * to reference into the cp_R, h_RT, and s_R arrays. + * + * @param temp Temperature (Kelvin) + * @param cp_R Vector of Dimensionless heat capacities. + * (length m_kk). + * @param h_RT Vector of Dimensionless enthalpies. + * (length m_kk). + * @param s_R Vector of Dimensionless entropies. + * (length m_kk). + */ + virtual void updatePropertiesTemp(const doublereal temp, + doublereal* cp_R, + doublereal* h_RT, doublereal* s_R) const ; - void updatePropertiesTemp(const doublereal temp, - doublereal* cp_R, - doublereal* h_RT, - doublereal* s_R) const ; + //!This utility function reports back the type of + //! parameterization and all of the parameters for the + //! species, index. + /*! + * All parameters are output variables + * + * @param n Species index + * @param type Integer type of the standard type + * @param tlow output - Minimum temperature + * @param thigh output - Maximum temperature + * @param pref output - reference pressure (Pa). + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state. + */ + virtual void reportParameters(int &n, int &type, + doublereal &tlow, doublereal &thigh, + doublereal &pref, + doublereal* const coeffs) const; - /** - * report all of the parameters that make up this - * interpolation. - */ - void reportParameters(int &n, int &type, - doublereal &tlow, doublereal &thigh, - doublereal &pref, - doublereal* const coeffs) const; - protected: - /** - * Number of intervals in the interpolating linear - * approximation. Number of points is one more than the - * number of intervals. - */ - int m_numIntervals; - /** - * Value of the enthalpy at T = 298.15. - * This value is tied to the Heat of formation of - * the species at 298.15. - */ - doublereal m_H298; - /** - * Points at which the standard state chemical potential - * are given. - */ - vector_fp m_t0_int; + //! Modify parameters for the standard state + /*! + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state. + */ + virtual void modifyParameters(doublereal* coeffs); - /* - * Mu0's are primary input data. They aren't strictly - * needed, but are kept here for convenience. - */ - vector_fp m_mu0_R_int; - vector_fp m_h0_R_int; - vector_fp m_s0_R_int; - vector_fp m_cp0_R_int; - doublereal m_lowT, m_highT, m_Pref; - int m_index; + protected: + + /** + * Number of intervals in the interpolating linear + * approximation. Number of points is one more than the + * number of intervals. + */ + int m_numIntervals; + + /** + * Value of the enthalpy at T = 298.15. + * This value is tied to the Heat of formation of + * the species at 298.15. + */ + doublereal m_H298; + + /** + * Points at which the standard state chemical potential + * are given. + */ + vector_fp m_t0_int; + + /** + * Mu0's are primary input data. They aren't strictly + * needed, but are kept here for convenience. + */ + vector_fp m_mu0_R_int; + + //! Dimensionless Enthalpies at the temperature points + vector_fp m_h0_R_int; + + //! Entropy at the points + vector_fp m_s0_R_int; + + //! Heat capacity at the points + vector_fp m_cp0_R_int; + //! Limiting low temperature + doublereal m_lowT; + //! Limiting high temperature + doublereal m_highT; - private: + //! Reference pressure + doublereal m_Pref; - }; + //! Species index + int m_index; + + private: - void installMu0ThermoFromXML(std::string speciesName, - SpeciesThermo& sp, int k, - const XML_Node* Mu0Node_ptr); + //! process the coefficients + /*! + * Mu0Poly(): + * + * In the constructor, we calculate and store the + * piecewise linear approximation to the thermodynamic + * functions. + * + * @param coeffs coefficients. These are defined as follows: + * + * coeffs[0] = number of points (integer) + * 1 = H298(J/kmol) + * 2 = T1 (Kelvin) + * 3 = mu1 (J/kmol) + * 4 = T2 (Kelvin) + * 5 = mu2 (J/kmol) + * 6 = T3 (Kelvin) + * 7 = mu3 (J/kmol) + * ........ + */ + void processCoeffs(const doublereal * coeffs); + + }; + + //! Install a Mu0 polynomial thermodynamic reference state + /*! + * Install a Mu0 polynomial thermodynamic reference state property + * parameterization for species k into a SpeciesThermo instance, + * getting the information from an XML database. + * + * @param speciesName Name of the species + * @param sp Owning SpeciesThermo object + * @param k Species index + * @param Mu0Node_ptr Pointer to the XML element containing the + * Mu0 information. + * + * @ingroup spthermo + */ + void installMu0ThermoFromXML(std::string speciesName, + SpeciesThermo& sp, int k, + const XML_Node* Mu0Node_ptr); } #endif - - - - diff --git a/Cantera/src/SpeciesThermoInterpType.h b/Cantera/src/SpeciesThermoInterpType.h index fc2397b31..f8104dc02 100644 --- a/Cantera/src/SpeciesThermoInterpType.h +++ b/Cantera/src/SpeciesThermoInterpType.h @@ -76,13 +76,8 @@ namespace Cantera { * computed property values should be written. This method updates only one value in * each array. * - * Temperature Polynomial: - * tt[0] = t; - * tt[1] = t*t; - * tt[2] = m_t[1]*t; - * tt[3] = m_t[2]*t; - * tt[4] = 1.0/t; - * tt[5] = std::log(t); + * The form and length of the Temperature Polynomial may vary depending on the + * parameterization. * * @param tempPoly vector of temperature polynomials * @param cp_R Vector of Dimensionless heat capacities. diff --git a/tools/doc/Cantera.cfg.in b/tools/doc/Cantera.cfg.in index 693251fd5..21fb76914 100755 --- a/tools/doc/Cantera.cfg.in +++ b/tools/doc/Cantera.cfg.in @@ -106,6 +106,7 @@ FILE_PATTERNS = Kinetics.h Kinetics.cpp \ NasaThermo.h NasaPoly1.h NasaPoly2.h \ ShomateThermo.h ShomatePoly.h SimpleThermo.h \ GeneralSpeciesThermo.h GeneralSpeciesThermo.cpp \ + ConstCpPoly.h ConstCpPoly.cpp Mu0Poly.h Mu0Poly.cpp \ utilities.h \ VPStandardStateTP.h VPStandardStateTP.cpp \ SingleSpeciesTP.h SingleSpeciesTP.cpp \