diff --git a/bin/exp3to2.sh b/bin/exp3to2.sh new file mode 100755 index 000000000..61866b41f --- /dev/null +++ b/bin/exp3to2.sh @@ -0,0 +1,20 @@ +#/bin/sh +# +# This sed script replaces 3 character exponents +# starting with 0 with 2 characters +# e-0xx -> e-xx +# e0xx -> exx +# E-0xx -> E-xx +# E0xx -> Exx +# where +# x is a digit +# +# It takes one argument, the file to be operated on. +# And, it writes to standard out. It may be used to do a +# replacement in place. +# +cp $1 .exp.txt +cat .exp.txt | sed 's/\([eE]-\)\(0\)\([0-9][0-9]\)/\1\3/g' | \ + sed 's/\([eE]\)\(0\)\([0-9][0-9]\)/\1\3/g' | \ + sed 's/\([eE]+\)\(0\)\([0-9][0-9]\)/\1\3/g' +rm .exp.txt diff --git a/include/cantera/thermo/GeneralSpeciesThermo.h b/include/cantera/thermo/GeneralSpeciesThermo.h index 2feab4881..68054ba53 100644 --- a/include/cantera/thermo/GeneralSpeciesThermo.h +++ b/include/cantera/thermo/GeneralSpeciesThermo.h @@ -13,26 +13,27 @@ #include "SpeciesThermoMgr.h" #include "NasaPoly1.h" #include "Nasa9Poly1.h" +#include "StatMech.h" #include "speciesThermoTypes.h" namespace Cantera { -//! A species thermodynamic property manager for a phase. -/*! - * This is a general manager that can handle a wide variety - * of species thermodynamic polynomials for individual species. - * It is slow, however, because it recomputes the functions of - * temperature needed for each species. What it does is to create - * a vector of SpeciesThermoInterpType objects. - * - * @ingroup mgrsrefcalc - */ -class GeneralSpeciesThermo : public SpeciesThermo -{ + //! A species thermodynamic property manager for a phase. + /*! + * This is a general manager that can handle a wide variety + * of species thermodynamic polynomials for individual species. + * It is slow, however, because it recomputes the functions of + * temperature needed for each species. What it does is to create + * a vector of SpeciesThermoInterpType objects. + * + * @ingroup mgrsrefcalc + */ + class GeneralSpeciesThermo : public SpeciesThermo + { -public: + public: //! Constructor GeneralSpeciesThermo(); @@ -213,7 +214,7 @@ public: #endif -private: + private: //! Provide the SpeciesthermoInterpType object /*! * provide access to the SpeciesThermoInterpType object. @@ -225,7 +226,7 @@ private: */ SpeciesThermoInterpType* provideSTIT(size_t k); -protected: + protected: /** * This is the main unknown in the object. It is @@ -260,7 +261,7 @@ protected: friend class VPSSMgr; -}; + }; } diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h index 315014876..62c5a5862 100644 --- a/include/cantera/thermo/IdealGasPhase.h +++ b/include/cantera/thermo/IdealGasPhase.h @@ -424,14 +424,64 @@ public: * property manager. * @see SpeciesThermo */ - virtual doublereal cp_mole() const; + virtual doublereal cp_mole() const; - /** - * Molar heat capacity at constant volume. Units: J/kmol/K. - * For an ideal gas mixture, - * \f[ \hat c_v = \hat c_p - \hat R. \f] - */ - virtual doublereal cv_mole() const; + /** + * Molar heat capacity at constant volume. Units: J/kmol/K. + * For an ideal gas mixture, + * \f[ \hat c_v = \hat c_p - \hat R. \f] + */ + virtual doublereal cv_mole() const; + + /** + * @returns species translational/rotational specific heat at + * constant volume. Inferred from the species gas + * constant and number of translational/rotational + * degrees of freedom. The translational/rotational + * modes are assumed to be fully populated, and are + * given by + * \f[ + * C^{tr}_{v,s} \equiv \frac{\partial e^{tr}_s}{\partial T} = \frac{5}{2} R_s + * \f] + * for diatomic molecules and + * \f[ + * C^{tr}_{v,s} \equiv \frac{\partial e^{tr}_s}{\partial T} = \frac{3}{2} R_s + * \f] + * for atoms. + */ + virtual doublereal cv_tr(doublereal ) const; + + /** + * @returns species translational specific heat at constant volume. + * Since the translational modes are assumed to be fully populated + * this is simply + * \f[ + * C^{trans}_{v,s} \equiv \frac{\partial e^{trans}_s}{\partial T} = \frac{3}{2} R_s + * \f] + */ + virtual doublereal cv_trans() const; + + /** + * @returns species rotational specific heat at constant volume. + * By convention, we lump the translational/rotational components + * \f[ + * C^{tr}_{v,s} \equiv C^{trans}_{v,s} + C^{rot}_{v,s} + * \f] + * so then + * \f[ + * C^{rot}_{v,s} \equiv C^{tr}_{v,s} - C^{trans}_{v,s} + * \f] + */ + virtual doublereal cv_rot(double atomicity) const; + + /** + * @returns species vibrational specific heat at + * constant volume. This is defined as + * \f[ + * C^{vib}_{v,s} \equiv \frac{\partial e^{vib}_s}{\partial T_V} = \frac{R_s \theta_{vs}^2 \exp\left(\theta_{vs}/T_V\right)}{\left[\left(\exp\left(\theta_{vs}/T_V\right)-1\right)T_V\right]^2} + * \f] + */ + virtual doublereal cv_vib(int k, doublereal T) const; //@} diff --git a/include/cantera/thermo/StatMech.h b/include/cantera/thermo/StatMech.h new file mode 100644 index 000000000..943a9a9b1 --- /dev/null +++ b/include/cantera/thermo/StatMech.h @@ -0,0 +1,220 @@ +/** + * @file StatMech.h + * Header for a single-species standard state object derived + * from + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#ifndef CT_STATMECH_H +#define CT_STATMECH_H + +/* + * $Revision: 279 $ + * $Date: 2009-12-05 13:08:43 -0600 (Sat, 05 Dec 2009) $ + */ + + + +#include "cantera/base/global.h" +#include "SpeciesThermoInterpType.h" +#include "SpeciesThermoMgr.h" +#include +#include + +namespace Cantera { + + //! + /*! + * @ingroup spthermo + */ + class StatMech : public SpeciesThermoInterpType { + + public: + + //! Empty constructor + StatMech(); + + + //! 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. + */ + StatMech(int n, doublereal tlow, doublereal thigh, doublereal pref, + const doublereal* coeffs, std::string my_name); + + //! copy constructor + /*! + * @param b object to be copied + */ + StatMech(const StatMech& b); + + //! assignment operator + /*! + * @param b object to be copied + */ + StatMech& operator=(const StatMech& b); + + //! Destructor + virtual ~StatMech(); + + //! 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; + + //! Returns an integer representing the species index + virtual size_t speciesIndex() const; + + //! Build a series of maps for the properties needed for species + int buildmap(); + + //! 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: + * tt[0] = t; + * tt[1] = t*t; + * tt[2] = t*t*t; + * tt[3] = t*t*t*t; + * tt[4] = 1.0/t; + * tt[5] = 1.0/(t*t); + * tt[6] = std::log(t); + * + * @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). + */ + virtual 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. + * + * Temperature Polynomial: + * tt[0] = t; + * tt[1] = t*t; + * tt[2] = t*t*t; + * tt[3] = t*t*t*t; + * tt[4] = 1.0/t; + * tt[5] = 1.0/(t*t); + * tt[6] = std::log(t); + * + * @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; + + //!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. There are + * 12 of them, designed to be compatible + * with the multiple temperature formulation. + * coeffs[0] is equal to one. + * coeffs[1] is min temperature + * coeffs[2] is max temperature + * coeffs[3+i] from i =0,9 are the coefficients themselves + */ + virtual void reportParameters(size_t& 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); + + protected: + //! lowest valid temperature + doublereal m_lowT; + //! highest valid temperature + doublereal m_highT; + //! standard-state pressure + doublereal m_Pref; + //! species index + int m_index; + //! array of polynomial coefficients + vector_fp m_coeff; + + std::string sp_name; + + //*generic species struct that contains everything we need here + // achtung: add doxygen markup here + // achtung: convert doubles to realdoubles + struct species + { + //Nominal T-R Degrees of freedom (cv = cfs*k*T) + doublereal cfs; + + // Mol. Wt. Molecular weight (kg/kmol) + doublereal mol_weight; + + // number of vibrational temperatures necessary + int nvib; + + // Theta_v Characteristic vibrational temperature(s) (K) + doublereal theta[5]; + }; + + std::map name_map; + + }; + +} +#endif + diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index c24b33c5f..9335cc46d 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -298,7 +298,16 @@ public: return err("cv_mole"); } - + /** + * @returns species vibrational specific heat at + * constant volume. + * + */ + /// Molar heat capacity at constant volume. Units: J/kmol/K. + virtual doublereal cv_vib(int, double) const { + return err("cv_vib"); + } + /** * @} * @name Mechanical Properties diff --git a/include/cantera/thermo/speciesThermoTypes.h b/include/cantera/thermo/speciesThermoTypes.h index ed45d6000..e390053a1 100644 --- a/include/cantera/thermo/speciesThermoTypes.h +++ b/include/cantera/thermo/speciesThermoTypes.h @@ -59,6 +59,10 @@ //! This is implemented in the class Nasa9PolyMultiTempRegion in Nasa9Poly1MultiTempRegion #define NASA9MULTITEMP 513 +//! Properties derived from theoretical considerations +//! This is implemented in the class statmech in StatMech.h +#define STAT 111 + //! Surface Adsorbate Model for a species on a surface. //! This is implemented in the class Adsorbate. #define ADSORBATE 1024 diff --git a/include/cantera/transport.h b/include/cantera/transport.h index 600afcedb..3a216886e 100644 --- a/include/cantera/transport.h +++ b/include/cantera/transport.h @@ -13,5 +13,6 @@ #include "transport/DustyGasTransport.h" #include "transport/MultiTransport.h" #include "transport/MixTransport.h" +#include "transport/PecosTransport.h" #include "transport/LiquidTransport.h" #endif diff --git a/include/cantera/transport/PecosTransport.h b/include/cantera/transport/PecosTransport.h new file mode 100644 index 000000000..c5f107cfa --- /dev/null +++ b/include/cantera/transport/PecosTransport.h @@ -0,0 +1,295 @@ +/** + * @file PecosTransport.h + * Header file defining class PecosTransport + */ + +/* $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_PECOSTRAN_H +#define CT_PECOSTRAN_H + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +// STL includes +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +// Cantera includes +#include "TransportBase.h" +#include "cantera/numerics/DenseMatrix.h" + +namespace Cantera { + + + class GasTransportParams; + + /** + * + * Class PecosTransport implements mixture-averaged transport + * properties for ideal gas mixtures. + * + */ + class PecosTransport : public Transport { + + public: + + virtual ~PecosTransport() {} + + virtual int model() const { return cPecosTransport; } + + //! Viscosity of the mixture + /*! + * + */ + virtual doublereal viscosity(); + + virtual void getSpeciesViscosities(doublereal* visc) + { update_T(); updateViscosity_T(); copy(m_visc.begin(), m_visc.end(), visc); } + + //! Return the thermal diffusion coefficients + /*! + * For this approximation, these are all zero. + */ + virtual void getThermalDiffCoeffs(doublereal* const dt); + + /*! returns the mixture thermal conductivity + * + * This is computed using the lumped model, + * \f[ + * k = k^{tr} + k^{ve} + * \f] + * where, + * \f[ + * k^{tr}= 5/2 \mu_s C_{v,s}^{trans} + \mu_s C_{v,s}^{rot} + * \f] + * and, + * \f[ + * k^{ve}= \mu_s C_{v,s}^{vib} + \mu_s C_{v,s}^{elec} + * \f] + * + */ + virtual doublereal thermalConductivity(); + + virtual void getBinaryDiffCoeffs(const int ld, doublereal* const d); + + + //! Mixture-averaged diffusion coefficients [m^2/s]. + /*! + * For the single species case or the pure fluid case + * the routine returns the self-diffusion coefficient. + * This is need to avoid a Nan result in the formula + * below. + */ + virtual void getMixDiffCoeffs(doublereal* const d); + + //! Returns the mixture-averaged diffusion coefficients [m^2/s]. + //! These are the coefficients for calculating the molar diffusive fluxes + //! from the species mole fraction gradients, computed according to + //! Eq. 12.176 in "Chemically Reacting Flow": + //! + //! \f[ D_{km}^* = \frac{1-X_k}{\sum_{j \ne k}^K X_j/\mathcal{D}_{kj}} \f] + //! + //! @param[out] d vector of mixture-averaged diffusion coefficients for + //! each species, length m_nsp. + void getMixDiffCoeffsMole(doublereal* const d); + + //! Returns the mixture-averaged diffusion coefficients [m^2/s]. + //! These are the coefficients for calculating the diffusive mass fluxes + //! from the species mass fraction gradients, computed according to + //! Eq. 12.178 in "Chemically Reacting Flow": + //! + //! \f[ \frac{1}{D_{km}} = \sum_{j \ne k}^K \frac{X_j}{\mathcal{D}_{kj}} + + //! \frac{X_k}{1-Y_k} \sum_{j \ne k}^K \frac{Y_j}{\mathcal{D}_{kj}} \f] + //! + //! @param[out] d vector of mixture-averaged diffusion coefficients for + //! each species, length m_nsp. + void getMixDiffCoeffsMass(doublereal* const d); + + virtual void getMobilities(doublereal* const mobil); + virtual void update_T(); + virtual void update_C(); + + //! Get the species diffusive mass fluxes wrt to + //! the mass averaged velocity, + //! given the gradients in mole fraction and temperature + /*! + * Units for the returned fluxes are kg m-2 s-1. + * + * @param ndim Number of dimensions in the flux expressions + * @param grad_T Gradient of the temperature + * (length = ndim) + * @param ldx Leading dimension of the grad_X array + * (usually equal to m_nsp but not always) + * @param grad_X Gradients of the mole fraction + * Flat vector with the m_nsp in the inner loop. + * length = ldx * ndim + * @param ldf Leading dimension of the fluxes array + * (usually equal to m_nsp but not always) + * @param fluxes Output of the diffusive mass fluxes + * Flat vector with the m_nsp in the inner loop. + * length = ldx * ndim + */ + virtual void getSpeciesFluxes(int ndim, + const doublereal* grad_T, + int ldx, + const doublereal* grad_X, + int ldf, doublereal* fluxes); + + //! Initialize the transport object + /*! + * + * Here we change all of the internal dimensions to be sufficient. + * We get the object ready to do property evaluations. + * + * @param tr Transport parameters for all of the species + * in the phase. + * + */ + virtual bool initGas( GasTransportParams& tr ); + + + /** + * + * Reads the transport table specified (currently defaults to internal file) + * + * Reads the user-specified transport table, appending new species + * data and/or replacing default species data. + * + */ + void read_blottner_transport_table (); + + friend class TransportFactory; + + /** + * Return a structure containing all of the pertinent parameters + * about a species that was used to construct the Transport + * properties in this object. + * + * @param k Species number to obtain the properties from. + */ + struct GasTransportData getGasTransportData(int); + + protected: + + /// default constructor + PecosTransport(); + + private: + + //! Calculate the pressure from the ideal gas law + doublereal pressure_ig() const { + return (m_thermo->molarDensity() * GasConstant * + m_thermo->temperature()); + } + + // mixture attributes + int m_nsp; + doublereal m_tmin, m_tmax; + vector_fp m_mw; + + // polynomial fits + vector m_visccoeffs; + vector m_condcoeffs; + vector m_diffcoeffs; + vector_fp m_polytempvec; + + // blottner fits + //int species = 20; + double a[500], b[500], c[500]; + + // property values + DenseMatrix m_bdiff; + vector_fp m_visc; + vector_fp m_sqvisc; + vector_fp m_cond; + + vector_fp m_molefracs; + + vector > m_poly; + vector m_astar_poly; + vector m_bstar_poly; + vector m_cstar_poly; + vector m_om22_poly; + DenseMatrix m_astar; + DenseMatrix m_bstar; + DenseMatrix m_cstar; + DenseMatrix m_om22; + + DenseMatrix m_phi; // viscosity weighting functions + DenseMatrix m_wratjk, m_wratkj1; + + vector_fp m_zrot; + vector_fp m_crot; + vector_fp m_cinternal; + vector_fp m_eps; + vector_fp m_alpha; + vector_fp m_dipoleDiag; + + doublereal m_temp, m_logt, m_kbt, m_t14, m_t32; + doublereal m_sqrt_kbt, m_sqrt_t; + + vector_fp m_sqrt_eps_k; + DenseMatrix m_log_eps_k; + vector_fp m_frot_298; + vector_fp m_rotrelax; + + doublereal m_lambda; + doublereal m_viscmix; + + // work space + vector_fp m_spwork; + + void updateThermal_T(); + void updateViscosity_T(); + void updateCond_T(); + void updateSpeciesViscosities(); + void updateDiff_T(); + void correctBinDiffCoeffs(); + bool m_viscmix_ok; + bool m_viscwt_ok; + bool m_spvisc_ok; + bool m_diffmix_ok; + bool m_bindiff_ok; + bool m_abc_ok; + bool m_spcond_ok; + bool m_condmix_ok; + + int m_mode; + + DenseMatrix m_epsilon; + DenseMatrix m_diam; + DenseMatrix incl; + bool m_debug; + + // specific heats + vector_fp cv_rot; + vector_fp cp_R; + vector_fp cv_int; + + }; +} +#endif diff --git a/include/cantera/transport/TransportBase.h b/include/cantera/transport/TransportBase.h index 6f6f51d11..6aa0bf272 100644 --- a/include/cantera/transport/TransportBase.h +++ b/include/cantera/transport/TransportBase.h @@ -53,6 +53,7 @@ const int cAqueousTransport = 750; const int cSimpleTransport = 770; const int cRadiativeTransport = 800; const int cWaterTransport = 721; +const int cPecosTransport = 900; //! \endcond // forward reference diff --git a/src/thermo/GeneralSpeciesThermo.cpp b/src/thermo/GeneralSpeciesThermo.cpp index 0e0fd2aba..4895dcd5c 100644 --- a/src/thermo/GeneralSpeciesThermo.cpp +++ b/src/thermo/GeneralSpeciesThermo.cpp @@ -22,106 +22,111 @@ namespace Cantera { -/* - * Constructors - */ -GeneralSpeciesThermo::GeneralSpeciesThermo() : + /* + * Constructors + */ + GeneralSpeciesThermo::GeneralSpeciesThermo() : SpeciesThermo(), m_tlow_max(0.0), m_thigh_min(1.0E30), m_p0(OneAtm), m_kk(0) -{ + { m_tlow_max = 0.0; m_thigh_min = 1.0E30; -} + } -GeneralSpeciesThermo:: -GeneralSpeciesThermo(const GeneralSpeciesThermo& b) : + GeneralSpeciesThermo:: + GeneralSpeciesThermo(const GeneralSpeciesThermo& b) : m_tlow_max(b.m_tlow_max), m_thigh_min(b.m_thigh_min), m_kk(b.m_kk) -{ + { m_sp.resize(m_kk, 0); for (size_t k = 0; k < m_kk; k++) { - SpeciesThermoInterpType* bk = b.m_sp[k]; - if (bk) { - m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType(); - } + SpeciesThermoInterpType* bk = b.m_sp[k]; + if (bk) { + m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType(); + } } -} + } -GeneralSpeciesThermo& -GeneralSpeciesThermo::operator=(const GeneralSpeciesThermo& b) -{ + GeneralSpeciesThermo& + GeneralSpeciesThermo::operator=(const GeneralSpeciesThermo& b) + { if (&b != this) { - m_tlow_max = b.m_tlow_max; - m_thigh_min = b.m_thigh_min; + m_tlow_max = b.m_tlow_max; + m_thigh_min = b.m_thigh_min; - for (size_t k = 0; k < m_kk; k++) { - SpeciesThermoInterpType* sp = m_sp[k]; - if (sp) { - delete sp; - m_sp[k] = 0; - } - } - m_kk = b.m_kk; - m_sp.resize(m_kk, 0); - for (size_t k = 0; k < m_kk; k++) { - SpeciesThermoInterpType* bk = b.m_sp[k]; - if (bk) { - m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType(); - } - } + for (size_t k = 0; k < m_kk; k++) { + SpeciesThermoInterpType* sp = m_sp[k]; + if (sp) { + delete sp; + m_sp[k] = 0; + } + } + m_kk = b.m_kk; + m_sp.resize(m_kk, 0); + for (size_t k = 0; k < m_kk; k++) { + SpeciesThermoInterpType* bk = b.m_sp[k]; + if (bk) { + m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType(); + } + } } return *this; -} + } -GeneralSpeciesThermo::~GeneralSpeciesThermo() -{ + GeneralSpeciesThermo::~GeneralSpeciesThermo() + { for (size_t k = 0; k < m_kk; k++) { - SpeciesThermoInterpType* sp = m_sp[k]; - if (sp) { - delete sp; - m_sp[k] = 0; - } + SpeciesThermoInterpType* sp = m_sp[k]; + if (sp) { + delete sp; + m_sp[k] = 0; + } } -} + } -SpeciesThermo* -GeneralSpeciesThermo::duplMyselfAsSpeciesThermo() const -{ + SpeciesThermo* + GeneralSpeciesThermo::duplMyselfAsSpeciesThermo() const + { GeneralSpeciesThermo* gsth = new GeneralSpeciesThermo(*this); return (SpeciesThermo*) gsth; -} + } -/* - * Install parameterization for a species. - * @param index Species index - * @param type parameterization type - * @param c coefficients. The meaning of these depends on - * the parameterization. - */ -void GeneralSpeciesThermo::install(std::string name, - size_t index, - int type, - const doublereal* c, - doublereal minTemp, - doublereal maxTemp, - doublereal refPressure) -{ + /* + * Install parameterization for a species. + * @param index Species index + * @param type parameterization type + * @param c coefficients. The meaning of these depends on + * the parameterization. + */ + void GeneralSpeciesThermo::install(std::string name, + size_t index, + int type, + const doublereal* c, + doublereal minTemp, + doublereal maxTemp, + doublereal refPressure) + { /* * Resize the arrays if necessary, filling the empty * slots with the zero pointer. */ + + if(minTemp <= 0.0) + { + throw CanteraError("Error in GeneralSpeciesThermo.cpp", + " Cannot take 0 tmin as input. \n\n"); + } + if (index >= m_kk) { - m_sp.resize(index+1, 0); - m_kk = index+1; + m_sp.resize(index+1, 0); + m_kk = index+1; } - //AssertThrow(m_sp[index] == 0, - // "Index position isn't null, duplication of assignment: " + int2str(index)); //int nfreq = 3; /* @@ -130,68 +135,74 @@ void GeneralSpeciesThermo::install(std::string name, switch (type) { case NASA1: - m_sp[index] = new NasaPoly1(index, minTemp, maxTemp, - refPressure, c); - break; + m_sp[index] = new NasaPoly1(index, minTemp, maxTemp, + refPressure, c); + break; case SHOMATE1: - m_sp[index] = new ShomatePoly(index, minTemp, maxTemp, - refPressure, c); - break; + m_sp[index] = new ShomatePoly(index, minTemp, maxTemp, + refPressure, c); + break; case CONSTANT_CP: case SIMPLE: - m_sp[index] = new ConstCpPoly(index, minTemp, maxTemp, - refPressure, c); - break; + m_sp[index] = new ConstCpPoly(index, minTemp, maxTemp, + refPressure, c); + break; case MU0_INTERP: - m_sp[index] = new Mu0Poly(index, minTemp, maxTemp, - refPressure, c); - break; + m_sp[index] = new Mu0Poly(index, minTemp, maxTemp, + refPressure, c); + break; case SHOMATE2: - m_sp[index] = new ShomatePoly2(index, minTemp, maxTemp, - refPressure, c); - break; + m_sp[index] = new ShomatePoly2(index, minTemp, maxTemp, + refPressure, c); + break; case NASA2: - m_sp[index] = new NasaPoly2(index, minTemp, maxTemp, - refPressure, c); - break; + m_sp[index] = new NasaPoly2(index, minTemp, maxTemp, + refPressure, c); + break; + + case STAT: + m_sp[index] = new StatMech(index, minTemp, maxTemp, + refPressure, c, name); + break; + case ADSORBATE: - m_sp[index] = new Adsorbate(index, minTemp, maxTemp, - refPressure, c); - break; + m_sp[index] = new Adsorbate(index, minTemp, maxTemp, + refPressure, c); + break; default: - throw UnknownSpeciesThermoModel( - "GeneralSpeciesThermo::install", - "unknown species type", int2str(type)); - break; + throw UnknownSpeciesThermoModel( + "GeneralSpeciesThermo::install", + "unknown species type", int2str(type)); + break; } if (!m_sp[index]) { - cout << "Null m_sp... index = " << index << endl; - cout << "type = " << type << endl; + cout << "Null m_sp... index = " << index << endl; + cout << "type = " << type << endl; } m_tlow_max = max(minTemp, m_tlow_max); m_thigh_min = min(maxTemp, m_thigh_min); -} + } -// Install a new species thermodynamic property -// parameterization for one species. -/* - * @param stit_ptr Pointer to the SpeciesThermoInterpType object - * This will set up the thermo for one species - */ -void GeneralSpeciesThermo::install_STIT(SpeciesThermoInterpType* stit_ptr) -{ + // Install a new species thermodynamic property + // parameterization for one species. + /* + * @param stit_ptr Pointer to the SpeciesThermoInterpType object + * This will set up the thermo for one species + */ + void GeneralSpeciesThermo::install_STIT(SpeciesThermoInterpType* stit_ptr) + { /* * Resize the arrays if necessary, filling the empty * slots with the zero pointer. */ if (!stit_ptr) { - throw CanteraError("GeneralSpeciesThermo::install_STIT", - "zero pointer"); + throw CanteraError("GeneralSpeciesThermo::install_STIT", + "zero pointer"); } size_t index = stit_ptr->speciesIndex(); if (index >= m_kk) { - m_sp.resize(index+1, 0); - m_kk = index+1; + m_sp.resize(index+1, 0); + m_kk = index+1; } AssertThrow(m_sp[index] == 0, "Index position isn't null, duplication of assignment: " + int2str(index)); @@ -208,161 +219,177 @@ void GeneralSpeciesThermo::install_STIT(SpeciesThermoInterpType* stit_ptr) m_tlow_max = max(minTemp, m_tlow_max); m_thigh_min = min(maxTemp, m_thigh_min); -} + } -void GeneralSpeciesThermo::installPDSShandler(size_t k, PDSS* PDSS_ptr, - VPSSMgr* vpssmgr_ptr) -{ + void GeneralSpeciesThermo::installPDSShandler(size_t k, PDSS* PDSS_ptr, + VPSSMgr* vpssmgr_ptr) + { STITbyPDSS* stit_ptr = new STITbyPDSS(k, vpssmgr_ptr, PDSS_ptr); install_STIT(stit_ptr); -} + } -/** - * Update the properties for one species. - */ -void GeneralSpeciesThermo:: -update_one(size_t k, doublereal t, doublereal* cp_R, - doublereal* h_RT, doublereal* s_R) const -{ + /** + * Update the properties for one species. + */ + void GeneralSpeciesThermo:: + update_one(size_t k, doublereal t, doublereal* cp_R, + doublereal* h_RT, doublereal* s_R) const + { SpeciesThermoInterpType* sp_ptr = m_sp[k]; if (sp_ptr) { - sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R); + sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R); } -} + } -/** - * Update the properties for all species. - */ -void GeneralSpeciesThermo:: -update(doublereal t, doublereal* cp_R, - doublereal* h_RT, doublereal* s_R) const -{ + /** + * Update the properties for all species. + */ + void GeneralSpeciesThermo:: + update(doublereal t, doublereal* cp_R, + doublereal* h_RT, doublereal* s_R) const + { vector::const_iterator _begin, _end; _begin = m_sp.begin(); _end = m_sp.end(); SpeciesThermoInterpType* sp_ptr = 0; for (; _begin != _end; ++_begin) { - sp_ptr = *(_begin); - if (sp_ptr) { - sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R); - } - // else { - // writelog("General::update: sp_ptr is NULL!\n"); - //} + sp_ptr = *(_begin); + if (sp_ptr) { + sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R); + } + // else { + // writelog("General::update: sp_ptr is NULL!\n"); + //} } -} + } -/** - * This utility function reports the type of parameterization - * used for the species, index. - */ -int GeneralSpeciesThermo::reportType(size_t index) const -{ + /** + * This utility function reports the type of parameterization + * used for the species, index. + */ + int GeneralSpeciesThermo::reportType(size_t index) const + { SpeciesThermoInterpType* sp = m_sp[index]; if (sp) { - return sp->reportType(); + return sp->reportType(); } return -1; -} + } -/** - * This utility function reports back the type of - * parameterization and all of the parameters for the - * species, index. - * For the NASA object, there are 15 coefficients. - */ -void GeneralSpeciesThermo:: -reportParams(size_t index, int& type, doublereal* const c, - doublereal& minTemp, doublereal& maxTemp, doublereal& refPressure) const -{ + /** + * This utility function reports back the type of + * parameterization and all of the parameters for the + * species, index. + * For the NASA object, there are 15 coefficients. + */ + void GeneralSpeciesThermo:: + reportParams(size_t index, int& type, doublereal* const c, + doublereal& minTemp, doublereal& maxTemp, doublereal& refPressure) const + { SpeciesThermoInterpType* sp = m_sp[index]; size_t n; if (sp) { - sp->reportParameters(n, type, minTemp, maxTemp, - refPressure, c); - if (n != index) { - throw CanteraError("GeneralSpeciesThermo::reportParams", - "Internal error encountered"); - } + sp->reportParameters(n, type, minTemp, maxTemp, + refPressure, c); + if (n != index) { + throw CanteraError("GeneralSpeciesThermo::reportParams", + "Internal error encountered"); + } } else { - type = -1; + type = -1; } -} + } -/** - * Return the lowest temperature at which the thermodynamic - * parameterization is valid. If no argument is supplied, the - * value is the one for which all species parameterizations - * are valid. Otherwise, if an integer argument is given, the - * value applies only to the species with that index. - */ -doublereal GeneralSpeciesThermo::minTemp(size_t k) const -{ + // //! Modify parameters for the standard state + // /*! + // * @param index Species index + // * @param c Vector of coefficients used to set the + // * parameters for the standard state. + // */ + // void GeneralSpeciesThermo:: + // modifyParams(size_t index, doublereal* c) + // { + // SpeciesThermoInterpType* sp = m_sp[index]; + // if (sp) { + // sp->modifyParameters(c); + // } + // } + + + /** + * Return the lowest temperature at which the thermodynamic + * parameterization is valid. If no argument is supplied, the + * value is the one for which all species parameterizations + * are valid. Otherwise, if an integer argument is given, the + * value applies only to the species with that index. + */ + doublereal GeneralSpeciesThermo::minTemp(size_t k) const + { if (k == npos) { - return m_tlow_max; + return m_tlow_max; } else { - SpeciesThermoInterpType* sp = m_sp[k]; - if (sp) { - return sp->minTemp(); - } + SpeciesThermoInterpType* sp = m_sp[k]; + if (sp) { + return sp->minTemp(); + } } return m_tlow_max; -} + } -doublereal GeneralSpeciesThermo::maxTemp(size_t k) const -{ + doublereal GeneralSpeciesThermo::maxTemp(size_t k) const + { if (k == npos) { - return m_thigh_min; + return m_thigh_min; } else { - SpeciesThermoInterpType* sp = m_sp[k]; - if (sp) { - return sp->maxTemp(); - } + SpeciesThermoInterpType* sp = m_sp[k]; + if (sp) { + return sp->maxTemp(); + } } return m_thigh_min; -} + } -doublereal GeneralSpeciesThermo::refPressure(size_t k) const -{ + doublereal GeneralSpeciesThermo::refPressure(size_t k) const + { if (k == npos) { - return m_p0; + return m_p0; } else { - SpeciesThermoInterpType* sp = m_sp[k]; - if (sp) { - return sp->refPressure(); - } + SpeciesThermoInterpType* sp = m_sp[k]; + if (sp) { + return sp->refPressure(); + } } return m_p0; -} + } -SpeciesThermoInterpType* GeneralSpeciesThermo::provideSTIT(size_t k) -{ + SpeciesThermoInterpType* GeneralSpeciesThermo::provideSTIT(size_t k) + { return (m_sp[k]); -} + } #ifdef H298MODIFY_CAPABILITY -doublereal GeneralSpeciesThermo::reportOneHf298(int k) const -{ + doublereal GeneralSpeciesThermo::reportOneHf298(int k) const + { SpeciesThermoInterpType* sp_ptr = m_sp[k]; doublereal h = -1.0; if (sp_ptr) { - h = sp_ptr->reportHf298(0); + h = sp_ptr->reportHf298(0); } return h; -} + } -void GeneralSpeciesThermo::modifyOneHf298(const int k, const doublereal Hf298New) -{ + void GeneralSpeciesThermo::modifyOneHf298(const int k, const doublereal Hf298New) + { SpeciesThermoInterpType* sp_ptr = m_sp[k]; if (sp_ptr) { - sp_ptr->modifyOneHf298(k, Hf298New); + sp_ptr->modifyOneHf298(k, Hf298New); } -} + } #endif diff --git a/src/thermo/IdealGasPhase.cpp b/src/thermo/IdealGasPhase.cpp index 564cfc817..5a8a7fba3 100644 --- a/src/thermo/IdealGasPhase.cpp +++ b/src/thermo/IdealGasPhase.cpp @@ -155,6 +155,98 @@ doublereal IdealGasPhase::cv_mole() const return cp_mole() - GasConstant; } + /** + * @returns species translational/rotational specific heat at + * constant volume. + * + * Either: $5/2 R_s$ or $3/2 R_s$ for molecules/atoms. + * + */ + doublereal IdealGasPhase::cv_tr (doublereal atomicity) const + { + // k is the species number + int dum = 0; + int type = 0; + doublereal c[12]; + doublereal minTemp; + doublereal maxTemp; + doublereal refPressure; + + m_spthermo->reportParams(dum,type,c,minTemp,maxTemp,refPressure); + + if(type != 111) + { + throw CanteraError("Error in IdealGasPhase.cpp", + "cv_tr only supported for StatMech!. \n\n"); + + } + + // see reportParameters for specific details + return c[3]; + } + + /** + * @returns species translational specific heat at constant volume. + */ + doublereal IdealGasPhase::cv_trans () const + { return 1.5*GasConstant; } + + /** + * @returns species rotational specific heat at constant volume. + * + */ + doublereal IdealGasPhase::cv_rot (double atom) const + { return std::max(cv_tr(atom) - cv_trans(), 0.); } + + /** + * @returns species vibrational specific heat at + * constant volume. + * + * C^{vib}_{v,s} = \frac{\partial e^{vib}_{v,s} }{\partial T} + * + * The species vibration energy ($e^{vib}_{v,s}$) is: + * + * 0: atom + * + * Diatomic: + * \f[ + * \frac{R_s \theta_{v,s}}{e^{\theta_{v,s}/T}-1} + * \f] + * + * General Molecules: + * \f[ + * \sum_i \frac{R_s \theta_{v,s,i}}{e^{\theta_{v,s,i}/T}-1} + * \f] + * + */ + doublereal IdealGasPhase::cv_vib (const int k, const doublereal T) const + { + + // k is the species number + int dum = 0; + int type = 0; + doublereal c[12]; + doublereal minTemp; + doublereal maxTemp; + doublereal refPressure; + + c[0] = temperature(); + + m_spthermo->reportParams(dum,type,c,minTemp,maxTemp,refPressure); + + // basic sanity check + if(type != 111) + { + throw CanteraError("Error in IdealGasPhase.cpp", + "cv_vib only supported for StatMech!. \n\n"); + + } + + // see reportParameters for specific details + return c[4]; + + } + // Mechanical Equation of State ---------------------------- // Chemical Potentials and Activities ---------------------- diff --git a/src/thermo/Makefile.am b/src/thermo/Makefile.am index a703953d1..8994b193d 100644 --- a/src/thermo/Makefile.am +++ b/src/thermo/Makefile.am @@ -20,33 +20,10 @@ cc_sources = ConstCpPoly.cpp ConstDensityThermo.cpp DebyeHuckel.cpp \ VPSSMgrFactory.cpp VPSSMgr_ConstVol.cpp VPSSMgr_General.cpp \ VPSSMgr_IdealGas.cpp VPSSMgr_Water_ConstVol.cpp \ VPSSMgr_Water_HKFT.cpp VPStandardStateTP.cpp WaterProps.cpp \ - WaterPropsIAPWS.cpp WaterPropsIAPWSphi.cpp WaterSSTP.cpp + WaterPropsIAPWS.cpp WaterPropsIAPWSphi.cpp WaterSSTP.cpp \ + StatMech.cpp - -#Elements.cpp Phase.cpp RedlichKisterVPSSTP.cpp \ - ThermoPhase.cpp IdealGasPhase.cpp ConstDensityThermo.cpp \ - SpeciesThermoFactory.cpp ConstCpPoly.cpp Nasa9Poly1.cpp \ - Nasa9PolyMultiTempRegion.cpp PDSS_Water.cpp PDSS_HKFT.cpp \ - Mu0Poly.cpp GeneralSpeciesThermo.cpp SurfPhase.cpp \ - ThermoFactory.cpp SpeciesThermoInterpType.cpp \ - VPSSMgr.cpp VPSSMgrFactory.cpp VPSSMgr_General.cpp \ - IdealSolnGasVPSS.cpp MolalityVPSSTP.cpp VPStandardStateTP.cpp \ - VPSSMgr_IdealGas.cpp VPSSMgr_ConstVol.cpp PDSS_ConstVol.cpp \ - PDSS_IdealGas.cpp PDSS_SSVol.cpp DebyeHuckel.cpp PDSS.cpp \ - WaterProps.cpp WaterPropsIAPWS.cpp WaterPropsIAPWSphi.cpp \ - VPSSMgr_Water_HKFT.cpp VPSSMgr_Water_ConstVol.cpp \ - PDSS_IonsFromNeutral.cpp IonsFromNeutralVPSSTP.cpp \ - GibbsExcessVPSSTP.cpp LatticePhase.cpp IdealMolalSoln.cpp \ - HMWSoln.cpp HMWSoln_input.cpp WaterSSTP.cpp \ - MetalSHEelectrons.cpp \ - IdealSolidSolnPhase.cpp LatticeSolidPhase.cpp \ - SingleSpeciesTP.cpp MineralEQ3.cpp \ - PseudoBinaryVPSSTP.cpp MargulesVPSSTP.cpp \ - StoichSubstanceSSTP.cpp PureFluidPhase.cpp \ - StoichSubstance.cpp -#PecosGasPhase.cpp - AM_CPPFLAGS = -I$(top_builddir)/include AM_CXXFLAGS = $(AM_CPPFLAGS) diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp index 0e063c7b7..40abbd7cb 100644 --- a/src/thermo/SpeciesThermoFactory.cpp +++ b/src/thermo/SpeciesThermoFactory.cpp @@ -17,6 +17,7 @@ using namespace std; #include "cantera/thermo/Mu0Poly.h" #include "Nasa9PolyMultiTempRegion.h" #include "cantera/thermo/Nasa9Poly1.h" +#include "cantera/thermo/StatMech.h" #include "cantera/thermo/AdsorbateThermo.h" #include "cantera/thermo/SpeciesThermoMgr.h" @@ -644,6 +645,53 @@ static void installNasa9ThermoFromXML(std::string speciesName, } } +/** + * Install a stat mech based property solver + * for species k into a SpeciesThermo instance. + */ +static void installStatMechThermoFromXML(std::string speciesName, + SpeciesThermo& sp, int k, + const std::vector& tp) +{ + const XML_Node * fptr = tp[0]; + int nRegTmp = tp.size(); + int nRegions = 0; + vector_fp cPoly; + StatMech *np_ptr = 0; + std::vector regionPtrs; + doublereal tmin, tmax, pref = OneAtm; + + // Loop over all of the possible temperature regions + for (int i = 0; i < nRegTmp; i++) { + fptr = tp[i]; + if (fptr) { + if (fptr->name() == "StatMech") { + if (fptr->hasChild("floatArray")) { + + tmin = fpValue((*fptr)["Tmin"]); + tmax = fpValue((*fptr)["Tmax"]); + if ((*fptr).hasAttrib("P0")) { + pref = fpValue((*fptr)["P0"]); + } + if ((*fptr).hasAttrib("Pref")) { + pref = fpValue((*fptr)["Pref"]); + } + + getFloatArray(fptr->child("floatArray"), cPoly, false); + if (cPoly.size() != 0) { + throw CanteraError("installStatMechThermoFromXML", + "Expected no coeff: this is not a polynomial representation"); + } + } + } + } + } + // set properties + tmin = 0.1; + vector_fp coeffs(1); + coeffs[0] = 0.0; + (&sp)->install(speciesName, k, STAT, &coeffs[0], tmin, tmax, pref); +} //! Install a Adsorbate polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance. /*! @@ -746,8 +794,11 @@ void SpeciesThermoFactory::installThermoForSpecies } else if (f->name() == "Mu0") { installMu0ThermoFromXML(speciesNode["name"], spthermo, k, f); } else if (f->name() == "NASA9") { - installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); + installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); } + else if (f->name() == "StatMech") { + installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp); + } else if (f->name() == "adsorbate") { installAdsorbateThermoFromXML(speciesNode["name"], spthermo, k, *f); } @@ -762,7 +813,11 @@ void SpeciesThermoFactory::installThermoForSpecies installNasaThermoFromXML(speciesNode["name"], spthermo, k, f0, f1); } else if (f0->name() == "Shomate" && f1->name() == "Shomate") { installShomateThermoFromXML(speciesNode["name"], spthermo, k, f0, f1); - } else if (f0->name() == "NASA9" && f1->name() == "NASA9") { + } + else if (f0->name() == "StatMech") { + installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp); + } + else if (f0->name() == "NASA9" && f1->name() == "NASA9") { installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); } else { throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"], @@ -772,13 +827,17 @@ void SpeciesThermoFactory::installThermoForSpecies const XML_Node* f0 = tp[0]; if (f0->name() == "NASA9") { installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); - } else { - throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"], - "multiple"); + } + else if (f0->name() == "StatMech") { + installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp); + } + else { + throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"], + "multiple"); } } else { - throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"], - "multiple"); + throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"], + "multiple"); } } } diff --git a/src/thermo/StatMech.cpp b/src/thermo/StatMech.cpp new file mode 100644 index 000000000..f680df1ce --- /dev/null +++ b/src/thermo/StatMech.cpp @@ -0,0 +1,811 @@ +/** + * @file StatMech.cpp + * \link Cantera::SpeciesThermoInterpType SpeciesThermoInterpType\endlink + */ + +/* $Author: hkmoffa $ + * $Revision: 279 $ + * $Date: 2009-12-05 13:08:43 -0600 (Sat, 05 Dec 2009) $ + */ +// Copyright 2007 Sandia National Laboratories + +#include "cantera/thermo/StatMech.h" +#include +#include + +namespace Cantera +{ + + // Statistical mechanics + /* + * @ingroup spthermo + */ + + //! Empty constructor + StatMech::StatMech() + : m_lowT(0.1), m_highT (1.0), + m_Pref(1.0E5), m_index (0) {} + + + // 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. + */ + StatMech::StatMech(int n, doublereal tlow, doublereal thigh, + doublereal pref, + const doublereal* coeffs, + std::string my_name) : + m_lowT (tlow), + m_highT (thigh), + m_Pref (pref), + m_index (n), + sp_name (my_name) + { + // should error on zero -- cannot take ln(0) + if(m_lowT <= 0.0){ + throw CanteraError("Error in StatMech.cpp", + " Cannot take 0 tmin as input. \n\n"); + } + buildmap(); + } + + // copy constructor + /* + * @param b object to be copied + */ + StatMech::StatMech(const StatMech& b) : + m_lowT (b.m_lowT), + m_highT (b.m_highT), + m_Pref (b.m_Pref), + m_index (b.m_index) + { + + } + + // assignment operator + /* + * @param b object to be copied + */ + StatMech& StatMech::operator=(const StatMech& b) { + if (&b != this) { + m_lowT = b.m_lowT; + m_highT = b.m_highT; + m_Pref = b.m_Pref; + m_index = b.m_index; + } + return *this; + } + + // Destructor + StatMech::~StatMech() { + } + + // duplicator + SpeciesThermoInterpType * + StatMech::duplMyselfAsSpeciesThermoInterpType() const { + StatMech* np = new StatMech(*this); + return (SpeciesThermoInterpType *) np; + } + + // Returns the minimum temperature that the thermo + // parameterization is valid + doublereal StatMech::minTemp() const + { + return m_lowT; + } + + // Returns the maximum temperature that the thermo + // parameterization is valid + doublereal StatMech::maxTemp() const { + return m_highT; + } + + // Returns the reference pressure (Pa) + doublereal StatMech::refPressure() const { return m_Pref; } + + // Returns an integer representing the type of parameterization + int StatMech::reportType() const { + return STAT; + } + + // Returns an integer representing the species index + size_t StatMech::speciesIndex() const { + return m_index; + } + + int StatMech::buildmap() + { + + // build vector of strings + std::vector SS; + + // now just iterate over name map to place each + // string in a key + + SS.push_back("Air"); + SS.push_back("CPAir"); + SS.push_back("Ar" ); + SS.push_back("Ar+" ); + SS.push_back("C" ); + SS.push_back("C+" ); + SS.push_back("C2" ); + SS.push_back("C2H" ); + SS.push_back("C2H2" ); + SS.push_back("C3" ); + SS.push_back("CF" ); + SS.push_back("CF2" ); + SS.push_back("CF3" ); + SS.push_back("CF4" ); + SS.push_back("CH" ); + SS.push_back("CH2" ); + SS.push_back("CH3" ); + SS.push_back("CH4" ); + SS.push_back("Cl" ); + SS.push_back("Cl2" ); + SS.push_back("CN" ); + SS.push_back("CN+" ); + SS.push_back("CO" ); + SS.push_back("CO+" ); + SS.push_back("CO2" ); + SS.push_back("F" ); + SS.push_back("F2" ); + SS.push_back("H" ); + SS.push_back("H+" ); + SS.push_back("H2" ); + SS.push_back("H2+" ); + SS.push_back("H2O" ); + SS.push_back("HCl" ); + SS.push_back("HCN" ); + SS.push_back("He" ); + SS.push_back("He+" ); + SS.push_back("N" ); + SS.push_back("N+" ); + SS.push_back("N2" ); + SS.push_back("CPN2" ); + SS.push_back("N2+" ); + SS.push_back("Ne" ); + SS.push_back("NCO" ); + SS.push_back("NH" ); + SS.push_back("NH+" ); + SS.push_back("NH2" ); + SS.push_back("NH3" ); + SS.push_back("NO" ); + SS.push_back("NO+" ); + SS.push_back("NO2" ); + SS.push_back("O" ); + SS.push_back("O+" ); + SS.push_back("O2" ); + SS.push_back("O2+" ); + SS.push_back("OH" ); + SS.push_back("Si" ); + SS.push_back("SiO" ); + SS.push_back("e"); + + // now place each species in a map + int ii; + for(ii=0; ii < SS.size(); ii++) + { + name_map[SS[ii]]=(new species); + + // init to crazy defaults + name_map[SS[ii]]->nvib = -1; + name_map[SS[ii]]->cfs = -1; + name_map[SS[ii]]->mol_weight = -1; + + name_map[SS[ii]]->theta[0] =0.0; + name_map[SS[ii]]->theta[1] =0.0; + name_map[SS[ii]]->theta[2] =0.0; + name_map[SS[ii]]->theta[3] =0.0; + name_map[SS[ii]]->theta[4] =0.0; + } + + // now set all species information + + // build Air + name_map["Air"]->cfs = 2.5; + name_map["Air"]->mol_weight=28.96; + name_map["Air"]->nvib=0; + + // build CPAir + name_map["CPAir"]->cfs = 2.5; + name_map["CPAir"]->mol_weight=28.96; + name_map["CPAir"]->nvib=0; + + // build Ar + name_map["Ar"]->cfs = 1.5; + name_map["Ar"]->mol_weight=39.944; + name_map["Ar"]->nvib=0; + + // build Ar+ + name_map["Ar+"]->cfs = 1.5; + name_map["Ar+"]->mol_weight=39.94345; + name_map["Ar+"]->nvib=0; + + // build C + name_map["C"]->cfs = 1.5; + name_map["C"]->mol_weight=12.011; + name_map["C"]->nvib=0; + + // build C+ + name_map["C+"]->cfs = 1.5; + name_map["C+"]->mol_weight=12.01045; + name_map["C+"]->nvib=0; + + // C2 + name_map["C2"]->cfs=2.5; + name_map["C2"]->mol_weight=24.022; + name_map["C2"]->nvib=1; + name_map["C2"]->theta[0]=2.6687e3; + + // C2H + name_map["C2H"]->cfs=2.5; + name_map["C2H"]->mol_weight=25.03; + name_map["C2H"]->nvib=3; + name_map["C2H"]->theta[0]=5.20100e+03; + name_map["C2H"]->theta[1]=7.20000e+03; + name_map["C2H"]->theta[2]=2.66100e+03; + + // C2H2 + name_map["C2H2"]->cfs=2.5; + name_map["C2H2"]->mol_weight=26.038; + name_map["C2H2"]->nvib=5; + name_map["C2H2"]->theta[0]=4.85290e+03; + name_map["C2H2"]->theta[1]=2.84000e+03; + name_map["C2H2"]->theta[2]=4.72490e+03; + name_map["C2H2"]->theta[3]=8.81830e+02; + name_map["C2H2"]->theta[4]=1.05080e+03; + + // C3 + name_map["C3"]->cfs=2.5; + name_map["C3"]->mol_weight=36.033; + name_map["C3"]->nvib=3; + name_map["C3"]->theta[0]=1.84500e+03; + name_map["C3"]->theta[1]=7.78700e+02; + name_map["C3"]->theta[2]=3.11760e+03; + + // CF + name_map["CF"]->cfs=2.5; + name_map["CF"]->mol_weight=31.00940; + name_map["CF"]->nvib=1; + name_map["CF"]->theta[0]=1.88214e+03; + + // CF2 + name_map["CF2"]->cfs=3; + name_map["CF2"]->mol_weight=50.00780; + name_map["CF2"]->nvib=3; + name_map["CF2"]->theta[0]=1.76120e+03; + name_map["CF2"]->theta[1]=9.56820e+02; + name_map["CF2"]->theta[2]=1.60000e+03; + + // CF3 + name_map["CF3"]->cfs=3; + name_map["CF3"]->mol_weight=69.00620; + name_map["CF3"]->nvib=4; + name_map["CF3"]->theta[0]=1.56800e+03; + name_map["CF3"]->theta[1]=1.00900e+03; + name_map["CF3"]->theta[2]=1.81150e+03; + name_map["CF3"]->theta[3]=7.36680e+02; + + // CF4 + name_map["CF4"]->cfs=3; + name_map["CF4"]->mol_weight=88.00460; + name_map["CF4"]->nvib=4; + name_map["CF4"]->theta[0]=1.30720e+03; + name_map["CF4"]->theta[1]=6.25892e+02; + name_map["CF4"]->theta[2]=1.84540e+03; + name_map["CF4"]->theta[3]=9.08950e+02; + + // CH + name_map["CH"]->cfs=2.5; + name_map["CH"]->mol_weight=13.01900; + name_map["CH"]->nvib=1; + name_map["CH"]->theta[0]=4.11290e+03; + + // CH2 + name_map["CH2"]->cfs=3; + name_map["CH2"]->mol_weight=14.02700; + name_map["CH2"]->nvib=3; + name_map["CH2"]->theta[0]=4.31650e+03; + name_map["CH2"]->theta[1]=1.95972e+03; + name_map["CH2"]->theta[2]=4.60432e+03; + + // CH3 + name_map["CH3"]->cfs=3; + name_map["CH3"]->mol_weight=15.03500; + name_map["CH3"]->nvib=4; + name_map["CH3"]->theta[0]=4.31650e+03; + name_map["CH3"]->theta[1]=8.73370e+02; + name_map["CH3"]->theta[2]=4.54960e+03; + name_map["CH3"]->theta[3]=2.01150e+03; + + // CH4 + name_map["CH4"]->cfs=3; + name_map["CH4"]->mol_weight=16.04300; + name_map["CH4"]->nvib=4; + name_map["CH4"]->theta[0]=4.19660e+03; + name_map["CH4"]->theta[1]=2.20620e+03; + name_map["CH4"]->theta[2]=4.34450e+03; + name_map["CH4"]->theta[3]=1.88600e+03; + + // Cl + name_map["Cl"]->cfs=1.5; + name_map["Cl"]->mol_weight=35.45300; + name_map["Cl"]->nvib=0; + + // Cl2 + name_map["Cl2"]->cfs=2.5; + name_map["Cl2"]->mol_weight=70.96; + name_map["Cl2"]->nvib=1; + name_map["Cl2"]->theta[0]=8.05355e+02; + + // CN + name_map["CN"]->cfs=2.5; + name_map["CN"]->mol_weight=26.01900; + name_map["CN"]->nvib=1; + name_map["CN"]->theta[0]=2.97610e+03; + + // CN+ + name_map["CN+"]->cfs=2.5; + name_map["CN+"]->mol_weight=26.01845; + name_map["CN+"]->nvib=1; + name_map["CN+"]->theta[0]=2.92520e+03; + + // CO + name_map["CO"]->cfs=2.5; + name_map["CO"]->mol_weight=28.01100; + name_map["CO"]->nvib=1; + name_map["CO"]->theta[0]=3.12200e+03; + + // CO+ + name_map["CO+"]->cfs=2.5; + name_map["CO+"]->mol_weight=28.01045; + name_map["CO+"]->nvib=1; + name_map["CO+"]->theta[0]=3.18800e+03; + + // CO2 + name_map["CO2"]->cfs=2.5; + name_map["CO2"]->mol_weight=44.01100; + name_map["CO2"]->nvib=3; + name_map["CO2"]->theta[0]=1.91870e+03; + name_map["CO2"]->theta[1]=9.59660e+02; + name_map["CO2"]->theta[2]=3.38210e+03; + + // F + name_map["F"]->cfs=1.5; + name_map["F"]->mol_weight=18.99840; + name_map["F"]->nvib=0; + + // F2 + name_map["F2"]->cfs=2.5; + name_map["F2"]->mol_weight=37.99680; + name_map["F2"]->nvib=1; + name_map["F2"]->theta[0]=1.32020e+03; + + // H + name_map["H"]->cfs=1.5; + name_map["H"]->mol_weight=1; + name_map["H"]->nvib=0; + + // H+ + name_map["H+"]->cfs=1.5; + name_map["H+"]->mol_weight=1.00745; + name_map["H+"]->nvib=0; + + // H2 + name_map["H2"]->cfs=2.5; + name_map["H2"]->mol_weight=2.01600; + name_map["H2"]->nvib=1; + name_map["H2"]->theta[0]=6.33140e+03; + + // H2+ + name_map["H2+"]->cfs=2.5; + name_map["H2+"]->mol_weight=2.01545; + name_map["H2+"]->nvib=1; + name_map["H2+"]->theta[0]=3.34280e+03; + + // H2O + name_map["H2O"]->cfs=3.0; + name_map["H2O"]->mol_weight=18.01600; + name_map["H2O"]->nvib=3; + name_map["H2O"]->theta[0]=5.26130e+03; + name_map["H2O"]->theta[1]=2.29460e+03; + name_map["H2O"]->theta[2]=5.40395e+03; + + // HCl + name_map["HCl"]->cfs=2.5; + name_map["HCl"]->mol_weight=36.46100; + name_map["HCl"]->nvib=1; + name_map["HCl"]->theta[0]=4.30330e+03; + + // HCN + name_map["HCN"]->cfs=2.5; + name_map["HCN"]->mol_weight=27.02700; + name_map["HCN"]->nvib=3; + name_map["HCN"]->theta[0]=3.01620e+03; + name_map["HCN"]->theta[1]=1.02660e+03; + name_map["HCN"]->theta[2]=4.76450e+03; + + // He + name_map["He"]->cfs=1.5; + name_map["He"]->mol_weight=4.00300; + name_map["He"]->nvib=0; + + // He+ + name_map["He+"]->cfs=1.5; + name_map["He+"]->mol_weight=4.00245; + name_map["He+"]->nvib=0; + + // N + name_map["N"]->cfs=1.5; + name_map["N"]->mol_weight=14.008; + name_map["N"]->nvib=0; + + // Ne + name_map["Ne"]->cfs=1.5; + name_map["Ne"]->mol_weight=20.17900; + name_map["Ne"]->nvib=0; + + // N+ + name_map["N+"]->cfs=1.5; + name_map["N+"]->mol_weight=14.00745; + name_map["N+"]->nvib=0; + + // N2 + name_map["N2"]->cfs=2.5; + name_map["N2"]->mol_weight=28.01600; + name_map["N2"]->nvib=1; + name_map["N2"]->theta[0]=3.39500e+03; + + // N2+ + name_map["N2+"]->cfs=2.5; + name_map["N2+"]->mol_weight=28.01545; + name_map["N2+"]->nvib=1; + name_map["N2+"]->theta[0]=3.17580e+03; + + // CPN2 + name_map["CPN2"]->cfs=2.5; + name_map["CPN2"]->mol_weight=28.01600; + name_map["CPN2"]->nvib=0; + + // NCO + name_map["NCO"]->cfs=2.5; + name_map["NCO"]->mol_weight=42.01900; + name_map["NCO"]->nvib=3; + name_map["NCO"]->theta[0]=1.83600e+03; + name_map["NCO"]->theta[1]=7.67100e+02; + name_map["NCO"]->theta[2]=2.76800e+03; + + // NH + name_map["NH"]->cfs=2.5; + name_map["NH"]->mol_weight=15.01600; + name_map["NH"]->nvib=1; + name_map["NH"]->theta[0]=4.72240e+03; + + // NH+ + name_map["NH+"]->cfs=2.5; + name_map["NH+"]->mol_weight=15.01545; + name_map["NH+"]->nvib=0; + + // NH2 + name_map["NH2"]->cfs=2.5; + name_map["NH2"]->mol_weight=16.02400; + name_map["NH2"]->nvib=0; + + // NH3 + name_map["NH3"]->cfs=2.5; + name_map["NH3"]->mol_weight=17.03200; + name_map["NH3"]->nvib=4; + name_map["NH3"]->theta[0]=4.78100e+03; + name_map["NH3"]->theta[1]=1.47040e+03; + name_map["NH3"]->theta[2]=4.95440e+03; + name_map["NH3"]->theta[3]=2.34070e+03; + + // NO + name_map["NO"]->cfs=2.5; + name_map["NO"]->mol_weight=30.00800; + name_map["NO"]->nvib=1; + name_map["NO"]->theta[0]=2.81700e+03; + + // NO+ + name_map["NO+"]->cfs=2.5; + name_map["NO+"]->mol_weight=30.00745; + name_map["NO+"]->nvib=1; + name_map["NO+"]->theta[0]=3.42100e+03; + + // NO2 + name_map["NO2"]->cfs=3; + name_map["NO2"]->mol_weight=46.00800; + name_map["NO2"]->nvib=3; + name_map["NO2"]->theta[0]=1.07900e+03; + name_map["NO2"]->theta[1]=1.90000e+03; + name_map["NO2"]->theta[2]=2.32700e+03; + + // O + name_map["O"]->cfs=1.5; + name_map["O"]->mol_weight=16.000; + name_map["O"]->nvib=0; + + // O+ + name_map["O+"]->cfs=1.5; + name_map["O+"]->mol_weight=15.99945; + name_map["O+"]->nvib=0; + + // O2 + name_map["O2"]->cfs=2.5; + name_map["O2"]->mol_weight=32.00000; + name_map["O2"]->nvib=1; + name_map["O2"]->theta[0]=2.23900e+03; + + // O2 + name_map["O2+"]->cfs=2.5; + name_map["O2+"]->mol_weight=31.99945; + name_map["O2+"]->nvib=1; + name_map["O2+"]->theta[0]=2.74120e+03; + + // OH + name_map["OH"]->cfs=2.5; + name_map["OH"]->mol_weight=17.00800; + name_map["OH"]->nvib=1; + name_map["OH"]->theta[0]=5.37820e+03; + + // Si + name_map["Si"]->cfs=1.5; + name_map["Si"]->mol_weight=28.08550; + name_map["Si"]->nvib=0; + + // SiO + name_map["SiO"]->cfs=2.5; + name_map["SiO"]->mol_weight=44.08550; + name_map["SiO"]->nvib=1; + name_map["SiO"]->theta[0]=1.78640e+03; + + // electron + name_map["e"]->cfs=1.5; + name_map["e"]->mol_weight=0.00055; + name_map["e"]->nvib=0; + + int dum = 0; + for(ii=0; ii < SS.size(); ii++) + { + // check nvib was initalized for all species + if(name_map[SS[ii]]->nvib == -1) + { + std::cout << name_map[SS[ii]]->nvib << std::endl; + throw CanteraError("Error in StatMech.cpp", + "nvib not initialized!. \n\n"); + + } + else + { + // check that theta is initalized + for(int i=0;invib;i++) + { + if(name_map[SS[ii]]->theta[i] <= 0.0) + { + throw CanteraError("Error in StatMech.cpp", + "theta not initalized!. \n\n"); + } + } + + // check that no non-zero theta exist + // for any theta larger than nvib! + for(int i=name_map[SS[ii]]->nvib;i<5;i++) + { + if(name_map[SS[ii]]->theta[i] != 0.0) + { + std::string err = "bad theta value for "+SS[ii]+"\n"; + throw CanteraError("StatMech.cpp",err); + } + } // done with for loop + } + + // check mol weight was initialized for all species + if(name_map[SS[ii]]->mol_weight == -1) + { + std::cout << name_map[SS[ii]]->mol_weight << std::endl; + throw CanteraError("Error in StatMech.cpp", + "mol_weight not initialized!. \n\n"); + + } + + // cfs was initialized for all species + if(name_map[SS[ii]]->cfs == -1) + { + std::cout << name_map[SS[ii]]->cfs << std::endl; + throw CanteraError("Error in StatMech.cpp", + "cfs not initialized!. \n\n"); + + } + + } // done with sanity checks + + // mark it zero, dude + return 0; + } + + // Update the properties for this species + /** + * + * \f[ + * \frac{C_p^0(T)}{R} = \frac{C_v^0(T)}{R} + 1 + * \f] + * + * Where, + * \f[ + * \frac{C_v^0(T)}{R} = \frac{C_v^{tr}(T)}{R} + \frac{C_v^{vib}(T)}{R} + * \f] + * + * + * @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 StatMech::updateProperties(const doublereal* tt, + doublereal* cp_R, doublereal* h_RT, + doublereal* s_R) const { + + std::map::iterator it; + + // get species name, to gather species properties + species* s; + + // pointer to map location of particular species + if(name_map.find(sp_name) != name_map.end()) + { + s = name_map.find(sp_name)->second; + } + else + { + //std::cout << sp_name << std::endl; + throw CanteraError("StatMech.cpp", + "species properties not found!. \n\n"); + } + + // translational + rotational specific heat + doublereal ctr = 0.0; + double theta = 0.0; + + // 5/2 * R for molecules, 3/2 * R for atoms + ctr += GasConstant * s->cfs; + + // vibrational energy + for(int i=0; i< s->nvib; i++) + { + theta = s->theta[i]; + ctr += GasConstant * theta * (theta* exp(theta/tt[0])/(tt[0]*tt[0]))/((exp(theta/tt[0])-1) * (exp(theta/tt[0])-1)); + } + + // Cp = Cv + R + doublereal cpdivR = ctr/GasConstant + 1; + + // ACTUNG: fix enthalpy and entropy + doublereal hdivRT = 0.0; + doublereal sdivR = 0.0; + + // return the computed properties in the location in the output + // arrays for this species + cp_R[m_index] = cpdivR; + h_RT[m_index] = hdivRT; + s_R [m_index] = sdivR; + } + + + // 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. + * + * Temperature Polynomial: + * tt[0] = t; + * tt[1] = t*t; + * tt[2] = t*t*t; + * tt[3] = t*t*t*t; + * tt[4] = 1.0/t; + * tt[5] = 1.0/(t*t); + * tt[6] = std::log(t); + * + * @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 StatMech::updatePropertiesTemp(const doublereal temp, + doublereal* cp_R, doublereal* h_RT, + doublereal* s_R) const { + double tPoly[1]; + tPoly[0] = temp; + updateProperties(tPoly, cp_R, h_RT, s_R); + } + + //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 species state data + */ + void StatMech::reportParameters(size_t &n, int &type, + doublereal &tlow, doublereal &thigh, + doublereal &pref, + doublereal* const coeffs) const + { + species* s; + + n = m_index; + type = STAT; + tlow = m_lowT; + thigh = m_highT; + pref = m_Pref; + for (int i = 0; i < 9; i++) + { + coeffs[i] = 0.0; + } + doublereal temp = coeffs[0]; + coeffs[1] = m_lowT; + coeffs[2] = m_highT; + + // get species name, to gather species properties + // pointer to map location of particular species + if(name_map.find(sp_name) != name_map.end()) + { + s = name_map.find(sp_name)->second; + } + else + { + //std::cout << sp_name << std::endl; + throw CanteraError("StatMech.cpp", + "species properties not found!. \n\n"); + } + + double theta = 0.0; + doublereal cvib = 0; + + // vibrational energy + for(int i=0; i< s->nvib; i++) + { + theta = s->theta[i]; + cvib += GasConstant * theta * (theta* exp(theta/temp)/(temp*temp))/((exp(theta/temp)-1) * (exp(theta/temp)-1)); + } + + // load vibrational energy + coeffs[3] = GasConstant * s->cfs; + coeffs[4] = cvib; + + } + + // Modify parameters for the standard state + /* + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state. + */ + void StatMech::modifyParameters(doublereal* coeffs) + { + + + + } + + +} + diff --git a/src/transport/Makefile.am b/src/transport/Makefile.am index 01bda90c9..614dbab6a 100644 --- a/src/transport/Makefile.am +++ b/src/transport/Makefile.am @@ -6,8 +6,8 @@ cc_sources = AqueousTransport.cpp LiquidTransport.cpp MMCollisionInt.cpp \ LiquidTranInteraction.cpp LiquidTransportData.cpp \ LiquidTransportParams.cpp TortuosityBase.cpp \ TortuosityBruggeman.cpp TortuosityMaxwell.cpp \ - TortuosityPercolation.cpp TransportParams.cpp GasTransport.cpp - + TortuosityPercolation.cpp TransportParams.cpp \ + GasTransport.cpp PecosTransport.cpp AM_CPPFLAGS = -I$(top_builddir)/include AM_CXXFLAGS = $(AM_CPPFLAGS) diff --git a/src/transport/PecosTransport.cpp b/src/transport/PecosTransport.cpp new file mode 100755 index 000000000..f6f342a51 --- /dev/null +++ b/src/transport/PecosTransport.cpp @@ -0,0 +1,758 @@ +/** + * @file PecosTransport.cpp + * Mixture-averaged transport properties. + * + */ + +/* $Author$ + * $Revision$ + * $Date$ + */ + +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/transport/PecosTransport.h" + +#include "cantera/base/utilities.h" +#include "cantera/transport/TransportParams.h" +#include "cantera/transport/TransportFactory.h" +#include "cantera/base/stringUtils.h" + +#include "cantera/thermo/IdealGasPhase.h" + +#include +using namespace std; + +/** + * Mole fractions below MIN_X will be set to MIN_X when computing + * transport properties. + */ +#define MIN_X 1.e-20 + +namespace Cantera { + + //////////////////// class PecosTransport methods ////////////// + + PecosTransport::PecosTransport() : + m_nsp(0), + m_tmin(-1.0), + m_tmax(100000.), + m_temp(-1.0), + m_logt(0.0) + { + + + } + + bool PecosTransport::initGas( GasTransportParams& tr ) { + + // constant substance attributes + m_thermo = tr.thermo; + m_nsp = m_thermo->nSpecies(); + m_tmin = m_thermo->minTemp(); + m_tmax = m_thermo->maxTemp(); + + // make a local copy of the molecular weights + m_mw.resize(m_nsp); + copy(m_thermo->molecularWeights().begin(), + m_thermo->molecularWeights().end(), m_mw.begin()); + + // copy polynomials and parameters into local storage + m_poly = tr.poly; + m_visccoeffs = tr.visccoeffs; + m_condcoeffs = tr.condcoeffs; + m_diffcoeffs = tr.diffcoeffs; + + m_zrot = tr.zrot; + m_crot = tr.crot; + m_epsilon = tr.epsilon; + m_mode = tr.mode_; + m_diam = tr.diam; + m_eps = tr.eps; + m_alpha = tr.alpha; + m_dipoleDiag.resize(m_nsp); + for (int i = 0; i < m_nsp; i++) { + m_dipoleDiag[i] = tr.dipole(i,i); + } + + m_phi.resize(m_nsp, m_nsp, 0.0); + m_wratjk.resize(m_nsp, m_nsp, 0.0); + m_wratkj1.resize(m_nsp, m_nsp, 0.0); + int j, k; + for (j = 0; j < m_nsp; j++) + for (k = j; k < m_nsp; k++) { + m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]); + m_wratjk(k,j) = sqrt(m_wratjk(j,k)); + m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]); + } + + m_polytempvec.resize(5); + m_visc.resize(m_nsp); + m_sqvisc.resize(m_nsp); + m_cond.resize(m_nsp); + m_bdiff.resize(m_nsp, m_nsp); + + m_molefracs.resize(m_nsp); + m_spwork.resize(m_nsp); + + // set flags all false + m_viscmix_ok = false; + m_viscwt_ok = false; + m_spvisc_ok = false; + m_spcond_ok = false; + m_condmix_ok = false; + m_spcond_ok = false; + m_diffmix_ok = false; + m_abc_ok = false; + + // read blottner fit parameters (A,B,C) + cout << "reading blottner"; + read_blottner_transport_table(); + cout << "done with blottner"; + + // set specific heats + cv_rot.resize(m_nsp); + cp_R.resize(m_nsp); + cv_int.resize(m_nsp); + + for (k = 0; k < m_nsp; k++) { + cv_rot[k] = tr.crot[k]; + cp_R[k] = ((IdealGasPhase*)tr.thermo)->cp_R_ref()[k]; + cv_int[k] = cp_R[k] - 2.5 - cv_rot[k]; + } + return true; + } + + + /********************************************************* + * + * Public methods + * + *********************************************************/ + + + /****************** viscosity ******************************/ + + /** + * The viscosity is computed using the Wilke mixture rule. + * \f[ + * \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}. + * \f] + * Here \f$ \mu_k \f$ is the viscosity of pure species \e k, + * and + * \f[ + * \Phi_{k,j} = \frac{\left[1 + * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2} + * {\sqrt{8}\sqrt{1 + M_k/M_j}} + * \f] + * @see updateViscosity_T(); + */ + doublereal PecosTransport::viscosity() { + + update_T(); + update_C(); + + if (m_viscmix_ok) return m_viscmix; + + doublereal vismix = 0.0; + int k; + // update m_visc and m_phi if necessary + if (!m_viscwt_ok) updateViscosity_T(); + + multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork)); + + for (k = 0; k < m_nsp; k++) { + vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom; + } + m_viscmix = vismix; + return vismix; + } + + /******************* binary diffusion coefficients **************/ + /* + * + * Using Ramshaw's self-consistent Effective Binary Diffusion + * (1990, J. Non-Equilib. Thermo) + * Adding more doxygen would be good here + */ + + void PecosTransport::getBinaryDiffCoeffs(const int ld, doublereal* const d) { + int i,j; + + update_T(); + + // if necessary, evaluate the binary diffusion coefficents + if (!m_bindiff_ok) updateDiff_T(); + + doublereal rp = 1.0/pressure_ig(); + for (i = 0; i < m_nsp; i++) + for (j = 0; j < m_nsp; j++) { + d[ld*j + i] = rp * m_bdiff(i,j); + } + } + + + void PecosTransport::getMobilities(doublereal* const mobil) { + int k; + getMixDiffCoeffs(DATA_PTR(m_spwork)); + doublereal c1 = ElectronCharge / (Boltzmann * m_temp); + for (k = 0; k < m_nsp; k++) { + mobil[k] = c1 * m_spwork[k] * m_thermo->charge(k); + } + } + + + /****************** thermal conductivity **********************/ + + /** + * The thermal conductivity is computed using the Wilke mixture rule. + * \f[ + * \k = \sum_s \frac{k_s X_s}{\sum_j \Phi_{s,j} X_j}. + * \f] + * Here \f$ \k_s \f$ is the conductivity of pure species \e s, + * and + * \f[ + * \Phi_{s,j} = \frac{\left[1 + * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_s}}\right)}\right]^2} + * {\sqrt{8}\sqrt{1 + M_s/M_j}} + * \f] + * @see updateCond_T(); + */ + doublereal PecosTransport::thermalConductivity() { + int k; + doublereal lambda = 0.0; + + update_T(); + update_C(); + + // update m_cond and m_phi if necessary + if (!m_spcond_ok) updateCond_T(); + if (!m_condmix_ok) { + + multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork)); + + for (k = 0; k < m_nsp; k++) { + lambda += m_molefracs[k] * m_cond[k]/m_spwork[k]; //denom; + } + + } + m_lambda = lambda; + return m_lambda; + + } + + + /****************** thermal diffusion coefficients ************/ + + /** + * Thermal diffusion is not considered in this pecos + * model. To include thermal diffusion, use transport manager + * MultiTransport instead. This methods fills out array dt with + * zeros. + */ + void PecosTransport::getThermalDiffCoeffs(doublereal* const dt) { + int k; + for (k = 0; k < m_nsp; k++) { + dt[k] = 0.0; + } + } + + /** + * @param ndim The number of spatial dimensions (1, 2, or 3). + * @param grad_T The temperature gradient (ignored in this model). + * @param ldx Leading dimension of the grad_X array. + * The diffusive mass flux of species \e k is computed from + * \f[ + * \vec{j}_k = -n M_k D_k \nabla X_k + \frac{\rho_k}{\rho} \sum_r n M_r D_r \nabla X_r + * \f] + * + * This is neglective pressure, forced and thermal diffusion. + * + */ + void PecosTransport::getSpeciesFluxes(int ndim, + const doublereal* grad_T, int ldx, const doublereal* grad_X, + int ldf, doublereal* fluxes) { + int n, k; + + update_T(); + update_C(); + + getMixDiffCoeffs(DATA_PTR(m_spwork)); + + const vector_fp& mw = m_thermo->molecularWeights(); + const doublereal* y = m_thermo->massFractions(); + doublereal rhon = m_thermo->molarDensity(); + + vector_fp sum(ndim,0.0); + + doublereal correction=0.0; + // grab 2nd (summation) term -- still need to multiply by mass fraction (\rho_s / \rho) + for (k = 0; k < m_nsp; k++) + { + correction += rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k]; + } + + for (n = 0; n < ndim; n++) { + for (k = 0; k < m_nsp; k++) { + fluxes[n*ldf + k] = -rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k] + y[k]*correction; + sum[n] += fluxes[n*ldf + k]; + } + } + // add correction flux to enforce sum to zero + for (n = 0; n < ndim; n++) { + for (k = 0; k < m_nsp; k++) { + fluxes[n*ldf + k] -= y[k]*sum[n]; + } + } + } + + /** + * Mixture-averaged diffusion coefficients [m^2/s]. + * + * For the single species case or the pure fluid case + * the routine returns the self-diffusion coefficient. + * This is need to avoid a Nan result in the formula + * below. + */ + void PecosTransport::getMixDiffCoeffs(doublereal* const d) { + + update_T(); + update_C(); + + // update the binary diffusion coefficients if necessary + if (!m_bindiff_ok) updateDiff_T(); + + int k, j; + doublereal mmw = m_thermo->meanMolecularWeight(); + doublereal sumxw = 0.0, sum2; + doublereal p = pressure_ig(); + if (m_nsp == 1) { + d[0] = m_bdiff(0,0) / p; + } else { + for (k = 0; k < m_nsp; k++) sumxw += m_molefracs[k] * m_mw[k]; + for (k = 0; k < m_nsp; k++) { + sum2 = 0.0; + for (j = 0; j < m_nsp; j++) { + if (j != k) { + sum2 += m_molefracs[j] / m_bdiff(j,k); + } + } + if (sum2 <= 0.0) { + d[k] = m_bdiff(k,k) / p; + } else { + d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2); + } + } + } + } + + void PecosTransport::getMixDiffCoeffsMole(doublereal* const d) + { + update_T(); + update_C(); + + // update the binary diffusion coefficients if necessary + if (!m_bindiff_ok) { + updateDiff_T(); + } + + doublereal p = m_thermo->pressure(); + if (m_nsp == 1) { + d[0] = m_bdiff(0,0) / p; + } else { + for (size_t k = 0; k < m_nsp; k++) { + double sum2 = 0.0; + for (size_t j = 0; j < m_nsp; j++) { + if (j != k) { + sum2 += m_molefracs[j] / m_bdiff(j,k); + } + } + if (sum2 <= 0.0) { + d[k] = m_bdiff(k,k) / p; + } else { + d[k] = (1 - m_molefracs[k]) / (p * sum2); + } + } + } + } + + void PecosTransport::getMixDiffCoeffsMass(doublereal* const d) + { + update_T(); + update_C(); + + // update the binary diffusion coefficients if necessary + if (!m_bindiff_ok) { + updateDiff_T(); + } + + doublereal mmw = m_thermo->meanMolecularWeight(); + doublereal p = m_thermo->pressure(); + + if (m_nsp == 1) { + d[0] = m_bdiff(0,0) / p; + } else { + for (size_t k=0; ktemperature(); + if (t == m_temp) return; + if (t <= 0.0) { + throw CanteraError("PecosTransport::update_T", + "negative temperature "+fp2str(t)); + } + m_temp = t; + m_logt = log(m_temp); + m_kbt = Boltzmann * m_temp; + m_sqrt_t = sqrt(m_temp); + m_t14 = sqrt(m_sqrt_t); + m_t32 = m_temp * m_sqrt_t; + m_sqrt_kbt = sqrt(Boltzmann*m_temp); + + // compute powers of log(T) + m_polytempvec[0] = 1.0; + m_polytempvec[1] = m_logt; + m_polytempvec[2] = m_logt*m_logt; + m_polytempvec[3] = m_logt*m_logt*m_logt; + m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt; + + // temperature has changed, so polynomial fits will need to be + // redone. + m_viscmix_ok = false; + m_spvisc_ok = false; + m_viscwt_ok = false; + m_spcond_ok = false; + m_diffmix_ok = false; + m_bindiff_ok = false; + m_abc_ok = false; + m_condmix_ok = false; + } + + /** + * @internal This is called the first time any transport property + * is requested from Mixture after the concentrations + * have changed. + */ + void PecosTransport::update_C() + { + // signal that concentration-dependent quantities will need to + // be recomputed before use, and update the local mole + // fractions. + + m_viscmix_ok = false; + m_diffmix_ok = false; + m_condmix_ok = false; + + m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); + + // add an offset to avoid a pure species condition + int k; + for (k = 0; k < m_nsp; k++) { + m_molefracs[k] = std::max(MIN_X, m_molefracs[k]); + } + } + + /************************************************************************* + * + * methods to update temperature-dependent properties + * + *************************************************************************/ + + /** + * + * Update the temperature-dependent parts of the mixture-averaged + * thermal conductivity. + * + * Calculated as, + * \f[ + * k= \mu_s (5/2 * C_{v,s}^{trans} + C_{v,s}^{rot} + C_{v,s}^{vib} + * \f] + * + * + */ + void PecosTransport::updateCond_T() { + + int k; + doublereal fivehalves = 5/2; + for (k = 0; k < m_nsp; k++) { + // need to add cv_elec in the future + m_cond[k] = m_visc[k] * ( fivehalves * cv_int[k] + cv_rot[k] + m_thermo->cv_vib(k,m_temp) ); + } + m_spcond_ok = true; + m_condmix_ok = false; + } + + + /** + * Update the binary diffusion coefficients. These are evaluated + * from the polynomial fits at unit pressure (1 Pa). + */ + void PecosTransport::updateDiff_T() { + + // evaluate binary diffusion coefficients at unit pressure + int i,j; + int ic = 0; + if (m_mode == CK_Mode) { + for (i = 0; i < m_nsp; i++) { + for (j = i; j < m_nsp; j++) { + m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic])); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } + else { + for (i = 0; i < m_nsp; i++) { + for (j = i; j < m_nsp; j++) { + m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec, + m_diffcoeffs[ic]); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } + + m_bindiff_ok = true; + m_diffmix_ok = false; + } + + + /** + * + * Update the pure-species viscosities. (Pa-s) = (kg/m/sec) + * + * Using Blottner fit for viscosity. Defines kinematic viscosity + * of the form + * \f[ + * \mu_s\left(T\right) = 0.10 \exp\left(A_s\left(\log T\right)^2 + B_s\log T + C_s\right) + * \f] + * where \f$ A_s \f$, \f$ B_s \f$, and \f$ C_s \f$ are constants. + * + */ + void PecosTransport::updateSpeciesViscosities() { + + // blottner + // return 0.10*std::exp(_a*(logT*logT) + _b*logT + _c); + + int k; + // iterate over species, update pure-species viscosity + for (k = 0; k < m_nsp; k++) { + m_visc[k] = 0.10*std::exp(a[k]*(m_logt*m_logt) + b[k]*m_logt + c[k]); + m_sqvisc[k] = sqrt(m_visc[k]); + } + + // time to update mixing + m_spvisc_ok = true; + } + + /* + * read_blottner_transport_table() + * loads up A B and C for blottner fits + * hardcoded for air, will need to generalize later + */ + + void PecosTransport::read_blottner_transport_table() + { + // istringstream blot + // ("Air 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + // "CPAir 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + // "N 1.15572000000e-02 6.03167900000e-01 -1.24327495000e+01\n" + // "N2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + // "CPN2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + // "NO 4.36378000000e-02 -3.35511000000e-02 -9.57674300000e+00\n" + // "O 2.03144000000e-02 4.29440400000e-01 -1.16031403000e+01\n" + // "O2 4.49290000000e-02 -8.26158000000e-02 -9.20194750000e+00\n" + // "C -8.3285e-3 0.7703240 -12.7378000\n" + // "C2 -8.4311e-3 0.7876060 -13.0268000\n" + // "C3 -8.4312e-3 0.7876090 -12.8240000\n" + // "C2H -2.4241e-2 1.0946550 -14.5835500\n" + // "CN -8.3811e-3 0.7860330 -12.9406000\n" + // "CO -0.019527394 1.013295 -13.97873\n" + // "CO2 -0.019527387 1.047818 -14.32212\n" + // "HCN -2.4241e-2 1.0946550 -14.5835500\n" + // "H -8.3912e-3 0.7743270 -13.6653000\n" + // "H2 -8.3346e-3 0.7815380 -13.5351000\n" + // "e 0.00000000000e+00 0.00000000000e+00 -1.16031403000e+01\n"); + + // + // from: AIAA-1997-2474 and Sandia Report SC-RR-70-754 + // + // # Air -- Identical to N2 fit + // # N -- Sandia Report SC-RR-70-754 + // # N2 -- Sandia Report SC-RR-70-754 + // # CPN2 -- Identical to N2 fit + // # NO -- Sandia Report SC-RR-70-754 + // # O -- Sandia Report SC-RR-70-754 + // # O2 -- Sandia Report SC-RR-70-754 + // # C -- AIAA-1997-2474 + // # C2 -- AIAA-1997-2474 + // # C3 -- AIAA-1997-2474 + // # C2H -- wild-ass guess: identical to HCN fit + // # CN -- AIAA-1997-2474 + // # CO -- AIAA-1997-2474 + // # CO2 -- AIAA-1997-2474 + // # HCN -- AIAA-1997-2474 + // # H -- AIAA-1997-2474 + // # H2 -- AIAA-1997-2474 + // # e -- Sandia Report SC-RR-70-754 + + istringstream blot + ("Air 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + "CPAir 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + "N 1.15572000000e-02 6.03167900000e-01 -1.24327495000e+01\n" + "N2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + "CPN2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" + "NO 4.36378000000e-02 -3.35511000000e-02 -9.57674300000e+00\n" + "O 2.03144000000e-02 4.29440400000e-01 -1.16031403000e+01\n" + "O2 4.49290000000e-02 -8.26158000000e-02 -9.20194750000e+00\n" + "C -8.3285e-3 0.7703240 -12.7378000\n" + "C2 -8.4311e-3 0.7876060 -13.0268000\n" + "C3 -8.4312e-3 0.7876090 -12.8240000\n" + "C2H -2.4241e-2 1.0946550 -14.5835500\n" + "CN -8.3811e-3 0.7860330 -12.9406000\n" + "CO -0.019527394 1.013295 -13.97873\n" + "CO2 -0.019527387 1.047818 -14.32212\n" + "HCN -2.4241e-2 1.0946550 -14.5835500\n" + "H -8.3912e-3 0.7743270 -13.6653000\n" + "H2 -8.3346e-3 0.7815380 -13.5351000\n" + "e 0.00000000000e+00 0.00000000000e+00 -1.16031403000e+01\n"); + + string line; + string name; + string ss1,ss2,ss3,ss4,sss; + int k; + int i = 0; + + while (std::getline(blot, line)) + { + + istringstream ss(line); + std::getline(ss, ss1, ' '); + std::getline(ss, ss2, ' '); + std::getline(ss, ss3, ' '); + std::getline(ss, ss4, ' '); + name = ss1; + + // now put coefficients in correct species + for (k = 0; k < m_nsp; k++) + { + string sss = m_thermo->speciesName(k); + + // this is the right species index + if(sss.compare(ss1) == 0) + { + a[k] = atof(ss2.c_str()); + b[k] = atof(ss3.c_str()); + c[k] = atof(ss4.c_str()); + + // index + i++; + } + else // default to air + { + + a[k] = 0.026; + b[k] = 0.3; + c[k] = -11.3; + } + + } // done with for loop + } + + + // for (k = 0; k < m_nsp; k++) + // { + // string sss = m_thermo->speciesName(k); + // cout << sss << endl; + // cout << a[k] << endl; + // cout << b[k] << endl; + // cout << c[k] << endl; + // } + + // simple sanity check + // if(i != m_nsp-1) + // { + // std::cout << "error\n" << i << std::endl; + // } + + } + + /** + * + * Update the temperature-dependent viscosity terms. + * Updates the array of pure species viscosities, and the + * weighting functions in the viscosity mixture rule. + * The flag m_visc_ok is set to true. + * + */ + void PecosTransport::updateViscosity_T() { + doublereal vratiokj, wratiojk, factor1; + + if (!m_spvisc_ok) updateSpeciesViscosities(); + + // see Eq. (9-5.15) of Reid, Prausnitz, and Poling + int j, k; + for (j = 0; j < m_nsp; j++) { + for (k = j; k < m_nsp; k++) { + vratiokj = m_visc[k]/m_visc[j]; + wratiojk = m_mw[j]/m_mw[k]; + + // Note that m_wratjk(k,j) holds the square root of + // m_wratjk(j,k)! + factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j); + m_phi(k,j) = factor1*factor1 / + (SqrtEight * m_wratkj1(j,k)); + m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk); + } + } + m_viscwt_ok = true; + } + + // /** + // * + // * This function returns a Transport data object for a given species. + // * + // */ + // struct GasTransportData PecosTransport:: + // getGasTransportData(int kSpecies) + // { + // struct GasTransportData td; + // td.speciesName = m_thermo->speciesName(kSpecies); + + // td.geometry = 2; + // if (m_crot[kSpecies] == 0.0) { + // td.geometry = 0; + // } else if (m_crot[kSpecies] == 1.0) { + // td.geometry = 1; + // } + // td.wellDepth = m_eps[kSpecies] / Boltzmann; + // td.dipoleMoment = m_dipoleDiag[kSpecies] * 1.0E25 / SqrtTen; + // td.diameter = m_diam(kSpecies, kSpecies) * 1.0E10; + // td.polarizability = m_alpha[kSpecies] * 1.0E30; + // td.rotRelaxNumber = m_zrot[kSpecies]; + + // return td; + // } + +} + diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index a26ad90c7..c1d01cd29 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -8,6 +8,7 @@ // known transport models #include "cantera/transport/MultiTransport.h" +#include "cantera/transport/PecosTransport.h" #include "cantera/transport/MixTransport.h" #include "cantera/transport/SolidTransport.h" #include "cantera/transport/DustyGasTransport.h" @@ -202,6 +203,7 @@ TransportFactory::TransportFactory() : m_models["Aqueous"] = cAqueousTransport; m_models["Simple"] = cSimpleTransport; m_models["User"] = cUserTransport; + m_models["Pecos"] = cPecosTransport; m_models["None"] = None; //m_models["Radiative"] = cRadiative; @@ -371,6 +373,11 @@ Transport* TransportFactory::newTransport(std::string transportModel, tr = new MixTransport; initTransport(tr, phase, CK_Mode, log_level); break; + // adding pecos transport model 2/13/12 + case cPecosTransport: + tr = new PecosTransport; + initTransport(tr, phase, 0, log_level); + break; case cSolidTransport: tr = new SolidTransport; tr->setThermo(*phase); diff --git a/test_problems/PecosTransport/PecosTransport.cpp b/test_problems/PecosTransport/PecosTransport.cpp index b4c97177a..25c3fa6d0 100644 --- a/test_problems/PecosTransport/PecosTransport.cpp +++ b/test_problems/PecosTransport/PecosTransport.cpp @@ -32,14 +32,18 @@ using namespace std; /*****************************************************************/ /*****************************************************************/ -#include "cantera/Cantera.h" #include "cantera/transport.h" #include "cantera/IdealGasMix.h" - #include "cantera/transport/TransportFactory.h" +//#include "Cantera.h" +//#include "transport.h" +//#include "IdealGasMix.h" + +//#include "TransportFactory.h" + using namespace Cantera; -using namespace Cantera_CXX; +//using namespace Cantera_CXX; void printDbl(double val) { if (fabs(val) < 5.0E-17) {