IdealGasPhase.cpp

Go to the documentation of this file.
00001 /**
00002  *
00003  *  @file IdealGasPhase.cpp
00004  *   ThermoPhase object for the ideal gas equation of
00005  * state - workhorse for %Cantera (see \ref thermoprops 
00006  * and class \link Cantera::IdealGasPhase IdealGasPhase\endlink).
00007  *
00008  */
00009 /*
00010  * $Id: IdealGasPhase.cpp 279 2009-12-05 19:08:43Z hkmoffa $
00011  */
00012 
00013 #ifdef WIN32
00014 #pragma warning(disable:4786)
00015 #pragma warning(disable:4503)
00016 #endif
00017 
00018 #include "ct_defs.h"
00019 #include "mix_defs.h"
00020 #include "IdealGasPhase.h"
00021 #include "SpeciesThermo.h"
00022 
00023 using namespace std;
00024 
00025 namespace Cantera {
00026   // Default empty Constructor
00027   IdealGasPhase::IdealGasPhase():
00028     m_mm(0),
00029     m_tmin(0.0),
00030     m_tmax(0.0),
00031     m_p0(-1.0),
00032     m_tlast(0.0),
00033     m_logc0(0.0)
00034   {
00035   }
00036 
00037   // Copy Constructor
00038   IdealGasPhase::IdealGasPhase(const IdealGasPhase& right):
00039     m_mm(right.m_mm),
00040     m_tmin(right.m_tmin),
00041     m_tmax(right.m_tmax),
00042     m_p0(right.m_p0),
00043     m_tlast(right.m_tlast),
00044     m_logc0(right.m_logc0)
00045   {
00046     /*
00047      * Use the assignment operator to do the brunt
00048      * of the work for the copy construtor.
00049      */
00050     *this = right;
00051   }
00052 
00053   // Asignment operator
00054   /*
00055    * Assignment operator for the object. Constructed
00056    * object will be a clone of this object, but will
00057    * also own all of its data.
00058    *
00059    * @param right Object to be copied.
00060    */
00061   IdealGasPhase& IdealGasPhase::
00062   operator=(const IdealGasPhase &right) {
00063     if (&right != this) {
00064       ThermoPhase::operator=(right);
00065       m_mm      = right.m_mm;
00066       m_tmin    = right.m_tmin;
00067       m_tmax    = right.m_tmax;
00068       m_p0      = right.m_p0;
00069       m_tlast   = right.m_tlast;
00070       m_logc0   = right.m_logc0;
00071       m_h0_RT   = right.m_h0_RT;
00072       m_cp0_R   = right.m_cp0_R;
00073       m_g0_RT   = right.m_g0_RT;
00074       m_s0_R    = right.m_s0_R;
00075       m_expg0_RT= right.m_expg0_RT;
00076       m_pe      = right.m_pe;
00077       m_pp      = right.m_pp;
00078     }
00079     return *this;
00080   }
00081 
00082   // Duplicator from the %ThermoPhase parent class
00083   /*
00084    * Given a pointer to a %ThermoPhase object, this function will
00085    * duplicate the %ThermoPhase object and all underlying structures.
00086    * This is basically a wrapper around the copy constructor.
00087    *
00088    * @return returns a pointer to a %ThermoPhase
00089    */
00090   ThermoPhase *IdealGasPhase::duplMyselfAsThermoPhase() const {
00091     ThermoPhase *igp = new IdealGasPhase(*this);
00092     return (ThermoPhase *) igp;
00093   }
00094 
00095   // Molar Thermodynamic Properties of the Solution ------------------
00096 
00097   /*
00098    * Molar internal energy. J/kmol. For an ideal gas mixture,
00099    * \f[
00100    * \hat u(T) = \sum_k X_k \hat h^0_k(T) - \hat R T,
00101    * \f]
00102    * and is a function only of temperature.
00103    * The reference-state pure-species enthalpies 
00104    * \f$ \hat h^0_k(T) \f$ are computed by the species thermodynamic 
00105    * property manager.
00106    * @see SpeciesThermo
00107    */
00108   doublereal IdealGasPhase::intEnergy_mole() const {
00109     return GasConstant * temperature()
00110       * ( mean_X(&enthalpy_RT_ref()[0]) - 1.0);
00111   }
00112 
00113   /*
00114    * Molar entropy. Units: J/kmol/K.
00115    * For an ideal gas mixture,
00116    * \f[
00117    * \hat s(T, P) = \sum_k X_k \hat s^0_k(T) - \hat R \log (P/P^0).
00118    * \f]
00119    * The reference-state pure-species entropies 
00120    * \f$ \hat s^0_k(T) \f$ are computed by the species thermodynamic 
00121    * property manager.
00122    * @see SpeciesThermo
00123    */
00124   doublereal IdealGasPhase::entropy_mole() const {
00125     return GasConstant * (mean_X(&entropy_R_ref()[0]) -
00126                           sum_xlogx() - std::log(pressure()/m_spthermo->refPressure()));
00127   }
00128 
00129   /*
00130    * Molar Gibbs free Energy for an ideal gas.
00131    * Units =  J/kmol.
00132    */
00133   doublereal IdealGasPhase::gibbs_mole() const {
00134     return enthalpy_mole() - temperature() * entropy_mole();
00135   }
00136 
00137   /*
00138    * Molar heat capacity at constant pressure. Units: J/kmol/K.
00139    * For an ideal gas mixture, 
00140    * \f[
00141    * \hat c_p(t) = \sum_k \hat c^0_{p,k}(T).
00142    * \f]
00143    * The reference-state pure-species heat capacities  
00144    * \f$ \hat c^0_{p,k}(T) \f$ are computed by the species thermodynamic 
00145    * property manager.
00146    * @see SpeciesThermo
00147    */
00148   doublereal IdealGasPhase::cp_mole() const {
00149     return GasConstant * mean_X(&cp_R_ref()[0]);
00150   }
00151 
00152   /*
00153    * Molar heat capacity at constant volume. Units: J/kmol/K.
00154    * For an ideal gas mixture,
00155    * \f[ \hat c_v = \hat c_p - \hat R. \f]
00156    */
00157   doublereal IdealGasPhase::cv_mole() const {
00158     return cp_mole() - GasConstant;
00159   }
00160 
00161   // Mechanical Equation of State ----------------------------
00162   // Chemical Potentials and Activities ----------------------
00163 
00164   /*
00165    * Returns the standard concentration \f$ C^0_k \f$, which is used to normalize
00166    * the generalized concentration.
00167    */
00168   doublereal IdealGasPhase::standardConcentration(int k) const {
00169     double p = pressure();
00170     return p/(GasConstant * temperature());
00171   }
00172 
00173   /*
00174    * Returns the natural logarithm of the standard 
00175    * concentration of the kth species
00176    */
00177   doublereal IdealGasPhase::logStandardConc(int k) const {
00178     _updateThermo();
00179     double p = pressure();
00180     double lc = std::log (p / (GasConstant * temperature()));
00181     return lc;
00182   }
00183 
00184   /* 
00185    * Get the array of non-dimensional activity coefficients 
00186    */
00187   void IdealGasPhase::getActivityCoefficients(doublereal *ac) const {
00188     for (int k = 0; k < m_kk; k++) {
00189       ac[k] = 1.0;
00190     }
00191   }
00192 
00193   /*
00194    * Get the array of chemical potentials at unit activity \f$
00195    * \mu^0_k(T,P) \f$.
00196    */
00197   void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const {
00198     const array_fp& gibbsrt = gibbs_RT_ref();
00199     scale(gibbsrt.begin(), gibbsrt.end(), muStar, _RT());
00200     double tmp = log (pressure() /m_spthermo->refPressure());
00201     tmp *=  GasConstant * temperature();
00202     for (int k = 0; k < m_kk; k++) {
00203       muStar[k] += tmp;  // add RT*ln(P/P_0)
00204     }
00205   }
00206 
00207   //  Partial Molar Properties of the Solution --------------
00208 
00209   void IdealGasPhase::getChemPotentials(doublereal* mu) const {
00210     getStandardChemPotentials(mu);
00211     //doublereal logp = log(pressure()/m_spthermo->refPressure());
00212     doublereal xx;
00213     doublereal rt = temperature() * GasConstant;
00214     //const array_fp& g_RT = gibbs_RT_ref();
00215     for (int k = 0; k < m_kk; k++) {
00216       xx = fmaxx(SmallNumber, moleFraction(k));
00217       mu[k] += rt*(log(xx));
00218     }
00219   }
00220   
00221   /*
00222    * Get the array of partial molar enthalpies of the species 
00223    * units = J / kmol
00224    */
00225   void IdealGasPhase::getPartialMolarEnthalpies(doublereal* hbar) const {
00226     const array_fp& _h = enthalpy_RT_ref();
00227     doublereal rt = GasConstant * temperature();
00228     scale(_h.begin(), _h.end(), hbar, rt);
00229   }
00230 
00231   /*
00232    * Get the array of partial molar entropies of the species 
00233    * units = J / kmol / K
00234    */
00235   void IdealGasPhase::getPartialMolarEntropies(doublereal* sbar) const {
00236     const array_fp& _s = entropy_R_ref();
00237     doublereal r = GasConstant;
00238     scale(_s.begin(), _s.end(), sbar, r);
00239     doublereal logp = log(pressure()/m_spthermo->refPressure());
00240     for (int k = 0; k < m_kk; k++) {
00241       doublereal xx = fmaxx(SmallNumber, moleFraction(k));
00242       sbar[k] += r * (- logp - log(xx));
00243     }
00244   }
00245 
00246   /*
00247    * Get the array of partial molar internal energies of the species 
00248    * units = J / kmol
00249    */
00250   void IdealGasPhase::getPartialMolarIntEnergies(doublereal* ubar) const {
00251     const array_fp& _h = enthalpy_RT_ref();
00252     doublereal rt = GasConstant * temperature();
00253     for (int k = 0; k < m_kk; k++) {
00254       ubar[k] =  rt * (_h[k] - 1.0);
00255     }
00256   }
00257 
00258   /*
00259    * Get the array of partial molar heat capacities
00260    */
00261   void IdealGasPhase::getPartialMolarCp(doublereal* cpbar) const {
00262     const array_fp& _cp = cp_R_ref();
00263     scale(_cp.begin(), _cp.end(), cpbar, GasConstant);
00264   }
00265 
00266   /*
00267    * Get the array of partial molar volumes
00268    * units = m^3 / kmol
00269    */
00270   void IdealGasPhase::getPartialMolarVolumes(doublereal* vbar) const {
00271     double vol = 1.0 / molarDensity();
00272     for (int k = 0; k < m_kk; k++) {
00273       vbar[k] = vol;
00274     }
00275   }
00276 
00277   // Properties of the Standard State of the Species in the Solution --
00278 
00279   /*
00280    * Get the nondimensional Enthalpy functions for the species
00281    * at their standard states at the current T and P of the
00282    * solution
00283    */
00284   void IdealGasPhase::getEnthalpy_RT(doublereal* hrt) const {
00285     const array_fp& _h = enthalpy_RT_ref();
00286     copy(_h.begin(), _h.end(), hrt);
00287   }
00288 
00289   /*
00290    * Get the array of nondimensional entropy functions for the 
00291    * standard state species
00292    * at the current <I>T</I> and <I>P</I> of the solution.
00293    */
00294   void IdealGasPhase::getEntropy_R(doublereal* sr) const {
00295     const array_fp& _s = entropy_R_ref();
00296     copy(_s.begin(), _s.end(), sr);
00297     double tmp = log (pressure() /m_spthermo->refPressure());
00298     for (int k = 0; k < m_kk; k++) {
00299       sr[k] -= tmp;
00300     }
00301   }
00302 
00303   /*
00304    * Get the nondimensional gibbs function for the species
00305    * standard states at the current T and P of the solution.
00306    */
00307   void IdealGasPhase::getGibbs_RT(doublereal* grt) const {
00308     const array_fp& gibbsrt = gibbs_RT_ref();
00309     copy(gibbsrt.begin(), gibbsrt.end(), grt);
00310     double tmp = log (pressure() /m_spthermo->refPressure());
00311     for (int k = 0; k < m_kk; k++) {
00312       grt[k] += tmp;
00313     }
00314   }
00315 
00316   /*
00317    * get the pure Gibbs free energies of each species assuming
00318    * it is in its standard state. This is the same as 
00319    * getStandardChemPotentials().
00320    */
00321   void IdealGasPhase::getPureGibbs(doublereal* gpure) const {
00322     const array_fp& gibbsrt = gibbs_RT_ref();
00323     scale(gibbsrt.begin(), gibbsrt.end(), gpure, _RT());
00324     double tmp = log (pressure() /m_spthermo->refPressure());
00325     tmp *= _RT();
00326     for (int k = 0; k < m_kk; k++) {
00327       gpure[k] += tmp;
00328     }
00329   }
00330 
00331   /*
00332    *  Returns the vector of nondimensional
00333    *  internal Energies of the standard state at the current temperature
00334    *  and pressure of the solution for each species.
00335    */
00336   void IdealGasPhase::getIntEnergy_RT(doublereal *urt) const {
00337     const array_fp& _h = enthalpy_RT_ref();
00338     for (int k = 0; k < m_kk; k++) {
00339       urt[k] = _h[k] - 1.0;
00340     }
00341   }
00342     
00343   /*
00344    * Get the nondimensional heat capacity at constant pressure
00345    * function for the species
00346    * standard states at the current T and P of the solution.
00347    */
00348   void IdealGasPhase::getCp_R(doublereal* cpr) const {
00349     const array_fp& _cpr = cp_R_ref();
00350     copy(_cpr.begin(), _cpr.end(), cpr);
00351   }
00352 
00353   /*
00354    *  Get the molar volumes of the species standard states at the current
00355    *  <I>T</I> and <I>P</I> of the solution.
00356    *  units = m^3 / kmol
00357    *
00358    * @param vol     Output vector containing the standard state volumes.
00359    *                Length: m_kk.
00360    */
00361   void IdealGasPhase::getStandardVolumes(doublereal *vol) const {
00362     double tmp = 1.0 / molarDensity();
00363     for (int k = 0; k < m_kk; k++) {
00364       vol[k] = tmp;
00365     }
00366   }
00367 
00368   // Thermodynamic Values for the Species Reference States ---------
00369 
00370   /*
00371    *  Returns the vector of nondimensional
00372    *  enthalpies of the reference state at the current temperature
00373    *  and reference presssure.
00374    */
00375   void IdealGasPhase::getEnthalpy_RT_ref(doublereal *hrt) const {
00376     const array_fp& _h = enthalpy_RT_ref();
00377     copy(_h.begin(), _h.end(), hrt);
00378   }
00379 
00380   /*
00381    *  Returns the vector of nondimensional
00382    *  enthalpies of the reference state at the current temperature
00383    *  and reference pressure.
00384    */
00385   void IdealGasPhase::getGibbs_RT_ref(doublereal *grt) const {
00386     const array_fp& gibbsrt = gibbs_RT_ref();
00387     copy(gibbsrt.begin(), gibbsrt.end(), grt);
00388   }
00389 
00390   /*
00391    *  Returns the vector of the 
00392    *  gibbs function of the reference state at the current temperature
00393    *  and reference pressure.
00394    *  units = J/kmol
00395    */
00396   void IdealGasPhase::getGibbs_ref(doublereal *g) const {
00397     const array_fp& gibbsrt = gibbs_RT_ref();
00398     scale(gibbsrt.begin(), gibbsrt.end(), g, _RT());
00399   }
00400 
00401   /*
00402    *  Returns the vector of nondimensional
00403    *  entropies of the reference state at the current temperature
00404    *  and reference pressure.
00405    */
00406   void IdealGasPhase::getEntropy_R_ref(doublereal *er) const {
00407     const array_fp& _s = entropy_R_ref();
00408     copy(_s.begin(), _s.end(), er);
00409   }
00410 
00411   /*
00412    *  Returns the vector of nondimensional
00413    *  internal Energies of the reference state at the current temperature
00414    *  of the solution and the reference pressure for each species.
00415    */
00416   void IdealGasPhase::getIntEnergy_RT_ref(doublereal *urt) const {
00417     const array_fp& _h = enthalpy_RT_ref();
00418     for (int k = 0; k < m_kk; k++) {
00419       urt[k] = _h[k] - 1.0;
00420     }
00421   }
00422 
00423   /*
00424    *  Returns the vector of nondimensional
00425    *  constant pressure heat capacities of the reference state
00426    *   at the current temperature and reference pressure.
00427    */
00428   void IdealGasPhase::getCp_R_ref(doublereal *cprt) const {
00429     const array_fp& _cpr = cp_R_ref();
00430     copy(_cpr.begin(), _cpr.end(), cprt);
00431   }
00432 
00433   void IdealGasPhase::getStandardVolumes_ref(doublereal *vol) const {
00434     doublereal tmp = _RT() / m_p0;
00435     for (int k = 0; k < m_kk; k++) {
00436       vol[k] = tmp;
00437     }
00438   }
00439 
00440 
00441     // new methods defined here -------------------------------
00442 
00443 
00444   void IdealGasPhase::initThermo() {
00445  
00446     m_mm = nElements();
00447     doublereal tmin = m_spthermo->minTemp();
00448     doublereal tmax = m_spthermo->maxTemp();
00449     if (tmin > 0.0) m_tmin = tmin;
00450     if (tmax > 0.0) m_tmax = tmax;
00451     m_p0 = refPressure();
00452 
00453     int leng = m_kk;
00454     m_h0_RT.resize(leng);
00455     m_g0_RT.resize(leng);
00456     m_expg0_RT.resize(leng);
00457     m_cp0_R.resize(leng);
00458     m_s0_R.resize(leng);
00459     m_pe.resize(leng, 0.0);
00460     m_pp.resize(leng);
00461   }
00462 
00463   /* 
00464    * Set mixture to an equilibrium state consistent with specified
00465    * chemical potentials and temperature. This method is needed by
00466    * the ChemEquil equillibrium solver.
00467    */
00468   void IdealGasPhase::setToEquilState(const doublereal* mu_RT) 
00469   {
00470     double tmp, tmp2;
00471     const array_fp& grt = gibbs_RT_ref();
00472 
00473     /*
00474      * Within the method, we protect against inf results if the
00475      * exponent is too high.
00476      *
00477      * If it is too low, we set
00478      * the partial pressure to zero. This capability is needed
00479      * by the elemental potential method.
00480      */
00481     doublereal pres = 0.0;
00482     for (int k = 0; k < m_kk; k++) {
00483       tmp = -grt[k] + mu_RT[k];
00484       if (tmp < -600.) {
00485         m_pp[k] = 0.0;
00486       } else if (tmp > 500.0) {
00487         tmp2 = tmp / 500.;
00488         tmp2 *= tmp2;
00489         m_pp[k] = m_p0 * exp(500.) * tmp2;
00490       } else {
00491         m_pp[k] = m_p0 * exp(tmp);
00492       }
00493       pres += m_pp[k];
00494     }
00495     // set state
00496     setState_PX(pres, &m_pp[0]);
00497   }
00498 
00499 
00500     /// This method is called each time a thermodynamic property is
00501     /// requested, to check whether the internal species properties
00502     /// within the object need to be updated.
00503     /// Currently, this updates the species thermo polynomial values
00504     /// for the current value of the temperature. A check is made
00505     /// to see if the temperature has changed since the last 
00506     /// evaluation. This object does not contain any persistent
00507     /// data that depends on the concentration, that needs to be
00508     /// updated. The state object modifies its concentration 
00509     /// dependent information at the time the setMoleFractions()
00510     /// (or equivalent) call is made.
00511     void IdealGasPhase::_updateThermo() const {
00512         doublereal tnow = temperature();
00513 
00514         // If the temperature has changed since the last time these
00515         // properties were computed, recompute them.
00516         if (m_tlast != tnow) {
00517             m_spthermo->update(tnow, &m_cp0_R[0], &m_h0_RT[0], 
00518                 &m_s0_R[0]);
00519             m_tlast = tnow;
00520 
00521             // update the species Gibbs functions
00522             int k;
00523             for (k = 0; k < m_kk; k++) {
00524                 m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
00525             }
00526             m_logc0 = log(m_p0/(GasConstant * tnow));
00527             m_tlast = tnow;
00528         }
00529     }
00530 }
00531 
Generated by  doxygen 1.6.3