From 39651964a78c6f2110bfbb3616412f12b9057204 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 15 Mar 2014 02:57:00 +0000 Subject: [PATCH] Added back Statmech and PecosTransport. Rarefied gas dynamics is an import application area for Sandia, and it represents a potentially important development direction for Cantera. These classes provide inroads into that area. --- include/cantera/thermo/GeneralSpeciesThermo.h | 1 + include/cantera/thermo/IdealGasPhase.h | 59 ++ include/cantera/thermo/StatMech.h | 171 +++++ include/cantera/thermo/speciesThermoTypes.h | 4 + include/cantera/transport.h | 1 + include/cantera/transport/PecosTransport.h | 314 +++++++++ include/cantera/transport/TransportBase.h | 1 + src/thermo/GeneralSpeciesThermo.cpp | 5 + src/thermo/IdealGasPhase.cpp | 55 ++ src/thermo/SpeciesThermoFactory.cpp | 53 ++ src/thermo/StatMech.cpp | 655 ++++++++++++++++++ src/transport/PecosTransport.cpp | 592 ++++++++++++++++ src/transport/TransportFactory.cpp | 8 + .../PecosTransport/PecosTransport.cpp | 252 +++++++ .../PecosTransport/output_blessed.txt | 176 +++++ test_problems/PecosTransport/runtest | 36 + test_problems/SConscript | 11 + test_problems/statmech/output_blessed.txt | 23 + test_problems/statmech/runtest_stat | 15 + test_problems/statmech/statmech_test_Fe.cpp | 57 ++ test_problems/statmech/test_stat_Fe.xml | 44 ++ .../statmech_properties/output_blessed.txt | 4 + .../statmech_properties.cpp | 61 ++ .../statmech_properties/test_stat.xml | 66 ++ .../statmech_test/output_blessed.txt | 0 test_problems/statmech_test/runtest_stat | 15 + test_problems/statmech_test/statmech_test.cpp | 95 +++ test_problems/statmech_test/test_stat.xml | 66 ++ .../statmech_test_poly/output_blessed.txt | 11 + .../statmech_test_poly/statmech_test_poly.cpp | 50 ++ .../statmech_test_poly/test_stat_err.xml | 47 ++ .../statmech_transport/output_blessed.txt | 1 + .../statmech_transport/statmech_transport.cpp | 103 +++ .../statmech_transport/test_stat_trans.xml | 100 +++ 34 files changed, 3152 insertions(+) create mode 100644 include/cantera/thermo/StatMech.h create mode 100644 include/cantera/transport/PecosTransport.h create mode 100644 src/thermo/StatMech.cpp create mode 100755 src/transport/PecosTransport.cpp create mode 100644 test_problems/PecosTransport/PecosTransport.cpp create mode 100644 test_problems/PecosTransport/output_blessed.txt create mode 100755 test_problems/PecosTransport/runtest create mode 100644 test_problems/statmech/output_blessed.txt create mode 100755 test_problems/statmech/runtest_stat create mode 100644 test_problems/statmech/statmech_test_Fe.cpp create mode 100644 test_problems/statmech/test_stat_Fe.xml create mode 100644 test_problems/statmech_properties/output_blessed.txt create mode 100644 test_problems/statmech_properties/statmech_properties.cpp create mode 100644 test_problems/statmech_properties/test_stat.xml create mode 100644 test_problems/statmech_test/output_blessed.txt create mode 100755 test_problems/statmech_test/runtest_stat create mode 100644 test_problems/statmech_test/statmech_test.cpp create mode 100644 test_problems/statmech_test/test_stat.xml create mode 100644 test_problems/statmech_test_poly/output_blessed.txt create mode 100644 test_problems/statmech_test_poly/statmech_test_poly.cpp create mode 100644 test_problems/statmech_test_poly/test_stat_err.xml create mode 100644 test_problems/statmech_transport/output_blessed.txt create mode 100644 test_problems/statmech_transport/statmech_transport.cpp create mode 100644 test_problems/statmech_transport/test_stat_trans.xml diff --git a/include/cantera/thermo/GeneralSpeciesThermo.h b/include/cantera/thermo/GeneralSpeciesThermo.h index b17a16aa3..8c945cfed 100644 --- a/include/cantera/thermo/GeneralSpeciesThermo.h +++ b/include/cantera/thermo/GeneralSpeciesThermo.h @@ -12,6 +12,7 @@ #include "SpeciesThermoMgr.h" #include "NasaPoly1.h" #include "Nasa9Poly1.h" +#include "StatMech.h" #include "speciesThermoTypes.h" namespace Cantera diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h index b216357da..8e4824d94 100644 --- a/include/cantera/thermo/IdealGasPhase.h +++ b/include/cantera/thermo/IdealGasPhase.h @@ -439,6 +439,65 @@ public: */ 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, + * \f[ + * C^{vib}_{v,s} = \frac{\partial e^{vib}_{v,s} }{\partial T} + * \f] + * where the species vibration energy \f$ e^{vib}_{v,s} \f$ is + * - atom: + * 0 + * - Diatomic: + * \f[ \frac{R_s \theta_{v,s}}{e^{\theta_{v,s}/T}-1} \f] + * - General Molecule: + * \f[ + * \sum_i \frac{R_s \theta_{v,s,i}}{e^{\theta_{v,s,i}/T}-1} + * \f] + */ + virtual doublereal cv_vib(int k, doublereal T) const; + //! @} //! @name Mechanical Equation of State //! @{ diff --git a/include/cantera/thermo/StatMech.h b/include/cantera/thermo/StatMech.h new file mode 100644 index 000000000..91f54bede --- /dev/null +++ b/include/cantera/thermo/StatMech.h @@ -0,0 +1,171 @@ +/** + * @file StatMech.h + * Header for a single-species standard state object derived + * from + */ +/* + * Copyright(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 + +#include "cantera/base/global.h" +#include "SpeciesThermoInterpType.h" +#include "SpeciesThermoMgr.h" + +namespace Cantera +{ + +//! Statistical mechanics +/*! + * @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, const 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); + + //! duplicator + virtual SpeciesThermoInterpType* + duplMyselfAsSpeciesThermoInterpType() const; + + //! Returns an integer representing the type of parameterization + virtual int reportType() 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. + * + * \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] + * + * Temperature Polynomial: + * tt[0] = 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. + * + * @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: + //! 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/speciesThermoTypes.h b/include/cantera/thermo/speciesThermoTypes.h index 4b5ab5d4d..9f7ebd88c 100644 --- a/include/cantera/thermo/speciesThermoTypes.h +++ b/include/cantera/thermo/speciesThermoTypes.h @@ -58,6 +58,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..e329b616c --- /dev/null +++ b/include/cantera/transport/PecosTransport.h @@ -0,0 +1,314 @@ +/** + * @file PecosTransport.h + * Header file defining class PecosTransport + */ + +// Copyright 2001 California Institute of Technology + +#ifndef CT_PECOSTRAN_H +#define CT_PECOSTRAN_H + +#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 int model() const { + return cPecosTransport; + } + + //! Viscosity of the mixture + /*! + * 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(); + */ + virtual doublereal viscosity(); + + virtual void getSpeciesViscosities(doublereal* const 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] + * + * 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(); + * @todo Reconcile these these formulas with the implementation + */ + virtual doublereal thermalConductivity(); + + //! binary diffusion coefficients + /*! + * Using Ramshaw's self-consistent Effective Binary Diffusion + * (1990, J. Non-Equilib. Thermo) + */ + virtual void getBinaryDiffCoeffs(const size_t 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. + */ + 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(); + + /** + * This is called the first time any transport property is requested from + * Mixture after the concentrations have changed. + */ + virtual void update_C(); + + //! Get the species diffusive mass fluxes wrt to the mass averaged + //! velocity, given the gradients in mole fraction and temperature + /*! + * 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 neglects pressure, forced and thermal diffusion. + * 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(size_t ndim, + const doublereal* const grad_T, + size_t ldx, + const doublereal* const grad_X, + size_t ldf, doublereal* const 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; + +protected: + 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; + vector_fp m_mw; + + // polynomial fits + std::vector m_visccoeffs; + std::vector m_condcoeffs; + std::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; + + std::vector > m_poly; + std::vector m_astar_poly; + std::vector m_bstar_poly; + std::vector m_cstar_poly; + std::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(); + + /** + * 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 updateViscosity_T(); + + /** + * 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 updateCond_T(); + + /** + * 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 updateSpeciesViscosities(); + + /** + * Update the binary diffusion coefficients. These are evaluated + * from the polynomial fits at unit pressure (1 Pa). + */ + 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 41f810949..77f22230b 100644 --- a/include/cantera/transport/TransportBase.h +++ b/include/cantera/transport/TransportBase.h @@ -54,6 +54,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 6a710fcaf..098cba653 100644 --- a/src/thermo/GeneralSpeciesThermo.cpp +++ b/src/thermo/GeneralSpeciesThermo.cpp @@ -145,6 +145,11 @@ void GeneralSpeciesThermo::install(const std::string& name, 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); diff --git a/src/thermo/IdealGasPhase.cpp b/src/thermo/IdealGasPhase.cpp index 6e5d193f1..7d69c0197 100644 --- a/src/thermo/IdealGasPhase.cpp +++ b/src/thermo/IdealGasPhase.cpp @@ -94,6 +94,61 @@ doublereal IdealGasPhase::cv_mole() const return cp_mole() - GasConstant; } +doublereal IdealGasPhase::cv_tr(doublereal atomicity) const +{ + // k is the species number + int dum = 0; + int type = m_spthermo->reportType(); + doublereal c[12]; + doublereal minTemp_; + doublereal maxTemp_; + doublereal refPressure_; + + if (type != 111) { + throw CanteraError("Error in IdealGasPhase.cpp", "cv_tr only supported for StatMech!. \n\n"); + } + + m_spthermo->reportParams(dum, type, c, minTemp_, maxTemp_, refPressure_); + + // see reportParameters for specific details + return c[3]; +} + +doublereal IdealGasPhase::cv_trans() const +{ + return 1.5 * GasConstant; +} + +doublereal IdealGasPhase::cv_rot(double atom) const +{ + return std::max(cv_tr(atom) - cv_trans(), 0.); +} + +doublereal IdealGasPhase::cv_vib(const int k, const doublereal T) const +{ + + // k is the species number + int dum = 0; + int type = m_spthermo->reportType(); + doublereal c[12]; + doublereal minTemp_; + doublereal maxTemp_; + doublereal refPressure_; + + c[0] = temperature(); + + // basic sanity check + if (type != 111) { + throw CanteraError("Error in IdealGasPhase.cpp", "cv_vib only supported for StatMech!. \n\n"); + } + + m_spthermo->reportParams(dum, type, c, minTemp_, maxTemp_, refPressure_); + + // see reportParameters for specific details + return c[4]; + +} + doublereal IdealGasPhase::standardConcentration(size_t k) const { double p = pressure(); diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp index 98d9e28cf..6e99dc3e6 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" @@ -578,6 +579,52 @@ static void installNasa9ThermoFromXML(const std::string& speciesName, } } +/** + * Install a stat mech based property solver + * for species k into a SpeciesThermo instance. + */ +static void installStatMechThermoFromXML(const std::string& speciesName, + SpeciesThermo& sp, int k, + const std::vector& tp) +{ + const XML_Node* fptr = tp[0]; + int nRegTmp = tp.size(); + vector_fp cPoly; + std::vector regionPtrs; + doublereal tmin = 0.0; + doublereal tmax = 0.0; + doublereal 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") { + tmin = fpValue((*fptr)["Tmin"]); + tmax = fpValue((*fptr)["Tmax"]); + if ((*fptr).hasAttrib("P0")) { + pref = fpValue((*fptr)["P0"]); + } + if ((*fptr).hasAttrib("Pref")) { + pref = fpValue((*fptr)["Pref"]); + } + if (fptr->hasChild("floatArray")) { + 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. /*! * This is called by method installThermoForSpecies if a Adsorbate block is found in the XML input. @@ -665,6 +712,8 @@ void SpeciesThermoFactory::installThermoForSpecies installMu0ThermoFromXML(speciesNode["name"], spthermo, k, f); } else if (f->name() == "NASA9") { 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); } else { @@ -678,6 +727,8 @@ 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() == "StatMech") { + installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp); } else if (f0->name() == "NASA9" && f1->name() == "NASA9") { installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); } else { @@ -688,6 +739,8 @@ void SpeciesThermoFactory::installThermoForSpecies const XML_Node* f0 = tp[0]; if (f0->name() == "NASA9") { installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); + } else if (f0->name() == "StatMech") { + installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp); } else { throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"], "multiple"); diff --git a/src/thermo/StatMech.cpp b/src/thermo/StatMech.cpp new file mode 100644 index 000000000..94213592e --- /dev/null +++ b/src/thermo/StatMech.cpp @@ -0,0 +1,655 @@ +/** + * @file StatMech.cpp + * \link Cantera::SpeciesThermoInterpType SpeciesThermoInterpType\endlink + */ + +// Copyright 2007 Sandia National Laboratories + +#include "cantera/thermo/StatMech.h" +#include + +namespace Cantera +{ +StatMech::StatMech() {} + +StatMech::StatMech(int n, doublereal tlow, doublereal thigh, + doublereal pref, + const doublereal* coeffs, + const std::string& my_name) : + SpeciesThermoInterpType(n, tlow, thigh, pref), + 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(); +} + +StatMech::StatMech(const StatMech& b) : + SpeciesThermoInterpType(b) +{ +} + +StatMech& StatMech::operator=(const StatMech& b) +{ + if (&b != this) { + SpeciesThermoInterpType::operator=(b); + } + return *this; +} + +SpeciesThermoInterpType* +StatMech::duplMyselfAsSpeciesThermoInterpType() const +{ + return new StatMech(*this); +} + +int StatMech::reportType() const +{ + return STAT; +} + +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 + size_t 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; + + for (ii=0; ii < SS.size(); ii++) { + // check nvib was initialized 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 initialized + for (int i=0; invib; i++) { + if (name_map[SS[ii]]->theta[i] <= 0.0) { + throw CanteraError("Error in StatMech.cpp", + "theta not initialized!. \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; +} + +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; +} + +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); +} + +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; + +} + +void StatMech::modifyParameters(doublereal* coeffs) +{ +} + +} diff --git a/src/transport/PecosTransport.cpp b/src/transport/PecosTransport.cpp new file mode 100755 index 000000000..5387a11e9 --- /dev/null +++ b/src/transport/PecosTransport.cpp @@ -0,0 +1,592 @@ +/** + * @file PecosTransport.cpp + * Mixture-averaged transport properties. + */ + +#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; + +namespace Cantera +{ + +PecosTransport::PecosTransport() : + m_nsp(0), + 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(); + + // 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) + read_blottner_transport_table(); + + // 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; +} + +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; +} + +void PecosTransport::getBinaryDiffCoeffs(const size_t 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); + } +} + +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; + +} + +void PecosTransport::getThermalDiffCoeffs(doublereal* const dt) +{ + int k; + for (k = 0; k < m_nsp; k++) { + dt[k] = 0.0; + } +} + +void PecosTransport::getSpeciesFluxes(size_t ndim, + const doublereal* const grad_T, + size_t ldx, const doublereal* const grad_X, + size_t ldf, doublereal* const fluxes) +{ + size_t n = 0; + int 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]; + } + } +} + +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 (int k = 0; k < m_nsp; k++) { + double sum2 = 0.0; + for (int 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 (int 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; +} + +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(Tiny, m_molefracs[k]); + } +} + +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; +} + +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; +} + +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; +} + +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] = fpValue(ss2); + b[k] = fpValue(ss3); + c[k] = fpValue(ss4); + + // 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; + // } + +} + +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; +} + +} diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index c03d2b75c..b632a7f2b 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" @@ -175,6 +176,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; for (map::iterator iter = m_models.begin(); @@ -341,7 +343,13 @@ Transport* TransportFactory::newTransport(const 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; initSolidTransport(tr, phase, log_level); tr->setThermo(*phase); diff --git a/test_problems/PecosTransport/PecosTransport.cpp b/test_problems/PecosTransport/PecosTransport.cpp new file mode 100644 index 000000000..ab60a6b45 --- /dev/null +++ b/test_problems/PecosTransport/PecosTransport.cpp @@ -0,0 +1,252 @@ +/** + * @file mixGasTransport.cpp + * test problem for mixture transport + */ + +// Example +// +// Test case for mixture transport in a gas +// The basic idea is to set up a gradient of some kind. +// Then the resulting transport coefficients out. +// Essentially all of the interface routines should be +// exercised and the results dumped out. +// +// A blessed solution test will make sure that the actual +// solution doesn't change as a function of time or +// further development. + +// perhaps, later, an analytical solution could be added + +#include "cantera/transport.h" +#include "cantera/IdealGasMix.h" +#include "cantera/transport/TransportFactory.h" + +#include + +using namespace std; +using namespace Cantera; + +void printDbl(double val) +{ + if (fabs(val) < 5.0E-17) { + cout << " nil"; + } else { + cout << val; + } +} + +int main(int argc, char** argv) +{ + size_t k; + string infile = "diamond.xml"; + + try { + + + IdealGasMix g("gri30.xml", "gri30_mix"); + size_t nsp = g.nSpecies(); + double pres = 1.0E5; + vector_fp Xset(nsp, 0.0); + Xset[0] = 0.269205 ; + Xset[1] = 0.000107082; + Xset[2] = 1.36377e-09 ; + Xset[3] = 4.35475e-10; + Xset[4] = 4.34036e-06 ; + Xset[5] = 0.192249; + Xset[6] = 3.59356e-13; + Xset[7] = 2.78061e-12 ; + Xset[8] = 4.7406e-18 ; + Xset[9] = 4.12955e-17 ; + Xset[10] = 2.58549e-14 ; + Xset[11] = 8.96502e-16 ; + Xset[12] = 6.09056e-11 ; + Xset[13] = 7.56752e-09 ; + Xset[14] = 0.192253; + Xset[15] = 0.0385036; + Xset[16] = 1.49596e-08 ; + Xset[17] = 2.22378e-08 ; + Xset[18] = 1.43096e-13 ; + Xset[19] = 1.45312e-15 ; + Xset[20] = 1.96948e-12 ; + Xset[21] = 8.41937e-19; + Xset[22] = 3.18852e-13 ; + Xset[23] = 7.93625e-18 ; + Xset[24] = 3.20653e-15 ; + Xset[25] = 1.15149e-19 ; + Xset[26] = 1.61189e-18 ; + Xset[27] = 1.4719e-15 ; + Xset[28] = 5.24728e-13 ; + Xset[29] = 6.90582e-17 ; + Xset[30] = 6.37248e-12 ; + Xset[31] =5.93728e-11 ; + Xset[32] = 2.71219e-09 ; + Xset[33] = 2.66645e-06 ; + Xset[34] = 6.57142e-11 ; + Xset[35] = 9.52453e-08 ; + Xset[36] = 1.26006e-14; + Xset[37] = 3.49802e-12; + Xset[38] = 1.19232e-11 ; + Xset[39] = 7.17782e-13 ; + Xset[40] = 1.85347e-07 ; + Xset[41] = 8.25325e-14 ; + Xset[42] = 5.00914e-20 ; + Xset[43] = 1.54407e-16 ; + Xset[44] =3.07176e-11 ; + Xset[45] =4.93198e-08 ; + Xset[46] =4.84792e-12 ; + Xset[47] = 0.307675 ; + Xset[48] =0; + Xset[49] =6.21649e-29; + Xset[50] = 8.42393e-28 ; + Xset[51] = 6.77865e-18; + Xset[52] = 2.19225e-16; + double T1 = 1500.; + + double sum = 0.0; + for (k = 0; k < nsp; k++) { + sum += Xset[k]; + } + for (k = 0; k < nsp; k++) { + Xset[k] /= sum; + } + + vector_fp X2set(nsp, 0.0); + X2set[0] = 0.25 ; + X2set[5] = 0.17; + X2set[14] = 0.15; + X2set[15] = 0.05; + X2set[47] = 0.38 ; + double T2 = 1200.; + + double dist = 0.1; + + vector_fp X3set(nsp, 0.0); + X3set[0] = 0.27 ; + X3set[5] = 0.15; + X3set[14] = 0.18; + X3set[15] = 0.06; + X3set[47] = 0.36 ; + double T3 = 1400.; + + vector_fp grad_T(3, 0.0); + + Array2D grad_X(nsp, 2, 0.0); + + + for (k = 0; k < nsp; k++) { + grad_X(k,0) = (X2set[k] - Xset[k])/dist; + grad_X(k,1) = (X3set[k] - Xset[k])/dist; + } + + grad_T[0] = (T2 - T1) / dist; + grad_T[1] = (T3 - T1) / dist; + + int log_level = 0; + Transport* tran = newTransportMgr("Pecos", &g, log_level=0); + PecosTransport* tranMix = dynamic_cast(tran); + + g.setState_TPX(1500.0, pres, DATA_PTR(Xset)); + + vector_fp mixDiffs(nsp, 0.0); + + tranMix->getMixDiffCoeffsMass(DATA_PTR(mixDiffs)); + printf(" Dump of the mixture Diffusivities:\n"); + for (k = 0; k < nsp; k++) { + string sss = g.speciesName(k); + printf(" %15s %13.5g\n", sss.c_str(), mixDiffs[k]); + } + + vector_fp specVisc(nsp, 0.0); + + tranMix->getSpeciesViscosities(DATA_PTR(specVisc)); + printf(" Dump of the species viscosities:\n"); + for (k = 0; k < nsp; k++) { + string sss = g.speciesName(k); + printf(" %15s %13.5g\n", sss.c_str(), specVisc[k]); + } + + vector_fp thermDiff(nsp, 0.0); + tranMix->getThermalDiffCoeffs(DATA_PTR(thermDiff)); + printf(" Dump of the Thermal Diffusivities :\n"); + for (k = 0; k < nsp; k++) { + string sss = g.speciesName(k); + printf(" %15s %13.5g\n", sss.c_str(), thermDiff[k]); + } + + printf("Viscoscity and thermal Cond vs. T\n"); + for (k = 0; k < 10; k++) { + T1 = 400. + 100. * k; + g.setState_TPX(T1, pres, DATA_PTR(Xset)); + double visc = tran->viscosity(); + double cond = tran->thermalConductivity(); + printf(" %13g %13.5g %13.5g\n", T1, visc, cond); + } + + g.setState_TPX(T1, pres, DATA_PTR(Xset)); + + Array2D Bdiff(nsp, nsp, 0.0); + printf("Binary Diffusion Coefficients H2 vs species\n"); + + tranMix->getBinaryDiffCoeffs(nsp, Bdiff.ptrColumn(0)); + for (k = 0; k < nsp; k++) { + string sss = g.speciesName(k); + printf(" H2 - %15s %13.5g %13.5g\n", sss.c_str(), Bdiff(0,k), Bdiff(k,0)); + } + + + vector_fp specMob(nsp, 0.0); + + tranMix->getMobilities(DATA_PTR(specMob)); + printf(" Dump of the species mobilities:\n"); + for (k = 0; k < nsp; k++) { + string sss = g.speciesName(k); + printf(" %15s %13.5g\n", sss.c_str(), specMob[k]); + } + + Array2D fluxes(nsp, 2, 0.0); + + tranMix->getSpeciesFluxes(2, DATA_PTR(grad_T), nsp, + grad_X.ptrColumn(0), nsp, fluxes.ptrColumn(0)); + printf(" Dump of the species fluxes:\n"); + double sum1 = 0.0; + double sum2 = 0.0; + double max1 = 0.0; + double max2 = 0.0; + for (k = 0; k < nsp; k++) { + string sss = g.speciesName(k); + printf(" %15s %13.5g %13.5g\n", sss.c_str(), fluxes(k,0), fluxes(k,1)); + sum1 += fluxes(k,0); + if (fabs(fluxes(k,0)) > max1) { + max1 = fabs(fluxes(k,0)); + } + sum2 += fluxes(k,1); + if (fabs(fluxes(k,1)) > max2) { + max2 = fabs(fluxes(k,0)); + } + } + + // Make sure roundoff error doesn't interfere with the printout. + // these should be zero. + if (fabs(sum1) * 1.0E14 > max1) { + printf("sum in x direction = %13.5g\n", sum1); + } else { + printf("sum in x direction = 0\n"); + } + if (fabs(sum2) * 1.0E14 > max2) { + printf("sum in y direction = %13.5g\n", sum1); + } else { + printf("sum in y direction = 0\n"); + } + + std::cout << "Sum of Diffusive Mass Fluxes: " << sum1 << std::endl; + std::cout << "Sum of Diffusive Mass Fluxes: " << sum2 << std::endl; + + + + } catch (CanteraError) { + showErrors(cout); + } + + return 0; +} +/***********************************************************/ diff --git a/test_problems/PecosTransport/output_blessed.txt b/test_problems/PecosTransport/output_blessed.txt new file mode 100644 index 000000000..4a360bb6a --- /dev/null +++ b/test_problems/PecosTransport/output_blessed.txt @@ -0,0 +1,176 @@ + Dump of the mixture Diffusivities: + H2 0.001204 + H 0.0022584 + O 0.00064032 + O2 0.00042264 + OH 0.00062946 + H2O 0.00051533 + HO2 0.00042019 + H2O2 0.00041763 + C 0.00059188 + CH 0.00068026 + CH2 0.0004608 + CH2(S) 0.0004608 + CH3 0.00045156 + CH4 0.00044985 + CO 0.0004216 + CO2 0.00034577 + HCO 0.00036192 + CH2O 0.00035926 + CH2OH 0.00035027 + CH3O 0.00035027 + CH3OH 0.00035149 + C2H 0.0003524 + C2H2 0.00034919 + C2H3 0.00034618 + C2H4 0.00034633 + C2H5 0.00031771 + C2H6 0.00031538 + HCCO 0.00053082 + CH2CO 0.00030634 + HCCOH 0.00030634 + N 0.00056491 + NH 0.00067404 + NH2 0.00066041 + NH3 0.00050222 + NNH 0.00040495 + NO 0.00041226 + NO2 0.00037026 + N2O 0.00033509 + HNO 0.00041903 + CN 0.00040648 + HCN 0.00035852 + H2CN 0.00035559 + HCNN 0.00053081 + HCNO 0.00033642 + HOCN 0.00033642 + HNCO 0.00033642 + NCO 0.00033783 + N2 0.00043046 + AR 0.00041248 + C3H7 0.00024676 + C3H8 0.00024577 + CH2CHO 0.00030507 + CH3CHO 0.00030384 + Dump of the species viscosities: + H2 4.4588e-05 + H 4.4588e-05 + O 4.4588e-05 + O2 4.4588e-05 + OH 4.4588e-05 + H2O 4.4588e-05 + HO2 4.4588e-05 + H2O2 4.4588e-05 + C 4.4588e-05 + CH 4.4588e-05 + CH2 4.4588e-05 + CH2(S) 4.4588e-05 + CH3 4.4588e-05 + CH4 4.4588e-05 + CO 4.4588e-05 + CO2 4.4588e-05 + HCO 4.4588e-05 + CH2O 4.4588e-05 + CH2OH 4.4588e-05 + CH3O 4.4588e-05 + CH3OH 4.4588e-05 + C2H 4.4588e-05 + C2H2 4.4588e-05 + C2H3 4.4588e-05 + C2H4 4.4588e-05 + C2H5 4.4588e-05 + C2H6 4.4588e-05 + HCCO 4.4588e-05 + CH2CO 4.4588e-05 + HCCOH 4.4588e-05 + N 4.4588e-05 + NH 4.4588e-05 + NH2 4.4588e-05 + NH3 4.4588e-05 + NNH 4.4588e-05 + NO 4.4588e-05 + NO2 4.4588e-05 + N2O 4.4588e-05 + HNO 4.4588e-05 + CN 4.4588e-05 + HCN 4.4588e-05 + H2CN 4.4588e-05 + HCNN 4.4588e-05 + HCNO 4.4588e-05 + HOCN 4.4588e-05 + HNCO 4.4588e-05 + NCO 4.4588e-05 + N2 4.4588e-05 + AR 4.4588e-05 + C3H7 4.4588e-05 + C3H8 4.4588e-05 + CH2CHO 4.4588e-05 + CH3CHO 4.4588e-05 + Dump of the Thermal Diffusivities : + H2 0 + H 0 + O 0 + O2 0 + OH 0 + H2O 0 + HO2 0 + H2O2 0 + C 0 + CH 0 + CH2 0 + CH2(S) 0 + CH3 0 + CH4 0 + CO 0 + CO2 0 + HCO 0 + CH2O 0 + CH2OH 0 + CH3O 0 + CH3OH 0 + C2H 0 + C2H2 0 + C2H3 0 + C2H4 0 + C2H5 0 + C2H6 0 + HCCO 0 + CH2CO 0 + HCCOH 0 + N 0 + NH 0 + NH2 0 + NH3 0 + NNH 0 + NO 0 + NO2 0 + N2O 0 + HNO 0 + CN 0 + HCN 0 + H2CN 0 + HCNN 0 + HCNO 0 + HOCN 0 + HNCO 0 + NCO 0 + N2 0 + AR 0 + C3H7 0 + C3H8 0 + CH2CHO 0 + CH3CHO 0 +Viscoscity and thermal Cond vs. T + + +************************************************ + Cantera Error! +************************************************ + + +Procedure: Error in IdealGasPhase.cpp +Error: cv_vib only supported for StatMech!. + + + + diff --git a/test_problems/PecosTransport/runtest b/test_problems/PecosTransport/runtest new file mode 100755 index 000000000..ab7e3d312 --- /dev/null +++ b/test_problems/PecosTransport/runtest @@ -0,0 +1,36 @@ +#!/bin/sh +# +# + +temp_success="1" +/bin/rm -f output.txt outputa.txt +tname="PecosTransport" +################################################################# +# +################################################################# +CANTERA_DATA=${CANTERA_DATA:=../../data/inputs}; export CANTERA_DATA + +CANTERA_BIN=${CANTERA_BIN:=../../bin} +./PecosTransport > output.txt +retnStat=$? +if [ $retnStat != "0" ] +then + temp_success="0" + echo "$tname ($tname test) returned with bad status, $retnStat, check output" + exit 1 +fi + +../../bin/exp3to2.sh output.txt > outputa.txt +diff -w outputa.txt output_blessed.txt > diff_test.out +retnStat=$? +if [ $retnStat = "0" ] +then + echo "successful diff comparison on $tname test" + exit 0 +else + echo "unsuccessful diff comparison on $tname test" + echo "FAILED" > csvCode.txt + temp_success="0" + exit 1 +fi + diff --git a/test_problems/SConscript b/test_problems/SConscript index 406bea101..f01cdd4aa 100644 --- a/test_problems/SConscript +++ b/test_problems/SConscript @@ -297,6 +297,7 @@ if haveConverters: Test('negA-cti', 'negATest', negA, 'negATest_blessed.out', options='noxNeg.cti', artifacts=negA_name) +CompileAndTest('pecosTransport', 'PecosTransport', 'pecosTransport', 'output_blessed.txt') CompileAndTest('printUtil', 'printUtilUnitTest', 'pUtest', 'output_blessed.txt') CompileAndTest('pureFluid', 'pureFluidTest', 'testPureWater', 'output_blessed.txt') if haveConverters: @@ -307,6 +308,16 @@ CompileAndTest('simpleTransport', 'simpleTransport', 'simpleTransport', 'output_blessed.txt') CompileAndTest('stoichSolidKinetics', 'stoichSolidKinetics', 'stoichSolidKinetics', 'output_blessed.txt') +CompileAndTest('statmech_test', 'statmech_test', + 'statmech_test', 'output_blessed.txt') +CompileAndTest('statmech_properties', 'statmech_properties', + 'statmech_properties', 'output_blessed.txt') +CompileAndTest('statmech_test_poly', 'statmech_test_poly', + 'statmech_test_poly', 'output_blessed.txt') +CompileAndTest('statmech_transport', 'statmech_transport', + 'statmech_transport', 'output_blessed.txt') +CompileAndTest('statmech_test_Fe', 'statmech', + 'statmech_test_Fe', 'output_blessed.txt') CompileAndTest('surfkin', 'surfkin', 'surfdemo', 'output_blessed.txt') CompileAndTest('surfSolver', 'surfSolverTest', 'surfaceSolver', None, arguments='haca2.xml', diff --git a/test_problems/statmech/output_blessed.txt b/test_problems/statmech/output_blessed.txt new file mode 100644 index 000000000..4d973a319 --- /dev/null +++ b/test_problems/statmech/output_blessed.txt @@ -0,0 +1,23 @@ + + +************************************************ + Cantera Error! +************************************************ + + +Procedure: StatMech.cpp +Error: species properties not found!. + + + +Procedure: StatMech.cpp +Error: species properties not found!. + + + +Procedure: StatMech.cpp +Error: species properties not found!. + + + + diff --git a/test_problems/statmech/runtest_stat b/test_problems/statmech/runtest_stat new file mode 100755 index 000000000..48da88d0c --- /dev/null +++ b/test_problems/statmech/runtest_stat @@ -0,0 +1,15 @@ +#!/bin/sh +# +# +temp_success="1" +/bin/rm -f output.txt outputa.txt +tname="mixGasTransport" +################################################################# +# +################################################################# +CANTERA_DATA=${CANTERA_DATA:=../../data/inputs}; export CANTERA_DATA + +CANTERA_BIN=${CANTERA_BIN:=../../bin} +./statmech_test > output.txt + +exit $? diff --git a/test_problems/statmech/statmech_test_Fe.cpp b/test_problems/statmech/statmech_test_Fe.cpp new file mode 100644 index 000000000..b89e4b7db --- /dev/null +++ b/test_problems/statmech/statmech_test_Fe.cpp @@ -0,0 +1,57 @@ +/** + * @file statmech + * test problem for statistical mechanics in cantera + */ + +// Example +// +// Test case to check error thrown if using Fe (not supported species) +// + +#include "cantera/transport.h" +#include "cantera/IdealGasMix.h" +#include "cantera/equil/equil.h" + +using namespace std; +using namespace Cantera; + +int main(int argc, char** argv) +{ + + try { + int k; + IdealGasMix g("test_stat_Fe.xml"); + int nsp = g.nSpecies(); + double pres = 1.0E5; + + vector_fp Xset(nsp, 0.0); + Xset[0] = 0.5 ; + Xset[1] = 0.5; + + g.setState_TPX(1500.0, pres, DATA_PTR(Xset)); + equilibrate(g, "TP", -1); + g.report(); + + vector_fp cp_R(nsp, 0.0); + g.getCp_R(DATA_PTR(cp_R)); + + for (int i=0; i + + + + + + + O H C N Na Cl Fe + + H Fe + + 2.165 + + + + + + + + + + + H:1 + + + + + 2.165 + + + + + Fe:1 + + + + + 2.165 + + + + + diff --git a/test_problems/statmech_properties/output_blessed.txt b/test_problems/statmech_properties/output_blessed.txt new file mode 100644 index 000000000..147c78b69 --- /dev/null +++ b/test_problems/statmech_properties/output_blessed.txt @@ -0,0 +1,4 @@ +2.5 +2.5 +2.5 +6.6558 diff --git a/test_problems/statmech_properties/statmech_properties.cpp b/test_problems/statmech_properties/statmech_properties.cpp new file mode 100644 index 000000000..df863dd22 --- /dev/null +++ b/test_problems/statmech_properties/statmech_properties.cpp @@ -0,0 +1,61 @@ +/** + * @file statmech + * test problem for statistical mechanics in cantera + */ + +// Example +// +// Test case for the statistical mechanics in cantera +// + +#include +#include +#include +#include +#include + +using namespace std; + +/*****************************************************************/ +/*****************************************************************/ + +#include "cantera/transport.h" +#include "cantera/IdealGasMix.h" +#include "cantera/equil/equil.h" + +#include "cantera/transport/TransportFactory.h" + +using namespace Cantera; + +int main(int argc, char** argv) +{ + + try { + int k; + IdealGasMix g("test_stat.xml"); + int nsp = g.nSpecies(); + double pres = 1.0E5; + + vector_fp Xset(nsp, 0.0); + Xset[0] = 0.5 ; + Xset[1] = 0.5; + + g.setState_TPX(1500.0, pres, DATA_PTR(Xset)); + equilibrate(g, "TP", -1); + + vector_fp cp_R(nsp, 0.0); + g.getCp_R(DATA_PTR(cp_R)); + + for (size_t i = 0; i < nsp; i++) { + cout << cp_R[i] << std::endl; + } + + } catch (CanteraError) { + showErrors(cout); + return 1; + } + + // Mark it zero! + return 0; + +} diff --git a/test_problems/statmech_properties/test_stat.xml b/test_problems/statmech_properties/test_stat.xml new file mode 100644 index 000000000..2938c4153 --- /dev/null +++ b/test_problems/statmech_properties/test_stat.xml @@ -0,0 +1,66 @@ + + + + + + + + + O H C N Na + + H O N NO2 + + 2.165 + + + + + + + + + + + H:1 + + + + + 2.165 + + + + + O:1 + + + + + 2.165 + + + + + N:1 + + + + + 2.165 + + + + + O:2 N:1 + + + + + + + 2.165 + + + + + diff --git a/test_problems/statmech_test/output_blessed.txt b/test_problems/statmech_test/output_blessed.txt new file mode 100644 index 000000000..e69de29bb diff --git a/test_problems/statmech_test/runtest_stat b/test_problems/statmech_test/runtest_stat new file mode 100755 index 000000000..48da88d0c --- /dev/null +++ b/test_problems/statmech_test/runtest_stat @@ -0,0 +1,15 @@ +#!/bin/sh +# +# +temp_success="1" +/bin/rm -f output.txt outputa.txt +tname="mixGasTransport" +################################################################# +# +################################################################# +CANTERA_DATA=${CANTERA_DATA:=../../data/inputs}; export CANTERA_DATA + +CANTERA_BIN=${CANTERA_BIN:=../../bin} +./statmech_test > output.txt + +exit $? diff --git a/test_problems/statmech_test/statmech_test.cpp b/test_problems/statmech_test/statmech_test.cpp new file mode 100644 index 000000000..039db7216 --- /dev/null +++ b/test_problems/statmech_test/statmech_test.cpp @@ -0,0 +1,95 @@ +/** + * @file statmech + * test problem for statistical mechanics in cantera + */ + +// Example +// +// Test case for the statistical mechanics in cantera +// + +#include "cantera/transport.h" +#include "cantera/IdealGasMix.h" +#include "cantera/equil/equil.h" + +using namespace std; +using namespace Cantera; + +int main(int argc, char** argv) +{ + + try { + int k; + IdealGasMix g("test_stat.xml"); + int nsp = g.nSpecies(); + double pres = 1.0E5; + + vector_fp Xset(nsp, 0.0); + Xset[0] = 0.5 ; + Xset[1] = 0.5; + + g.setState_TPX(1500.0, pres, DATA_PTR(Xset)); + equilibrate(g, "TP", -1); + + vector_fp cp_R(nsp, 0.0); + g.getCp_R(DATA_PTR(cp_R)); + + //for(int i=0;i= tol) { + double diff = cp_R[3]-sol; + std::cout << "Error for Species NO2!\n"; + std::cout << "Diff was: " << diff << "\n"; + return 1; + } + + } catch (CanteraError) { + showErrors(cout); + return 1; + } + + // Mark it zero! + return 0; + +} diff --git a/test_problems/statmech_test/test_stat.xml b/test_problems/statmech_test/test_stat.xml new file mode 100644 index 000000000..2938c4153 --- /dev/null +++ b/test_problems/statmech_test/test_stat.xml @@ -0,0 +1,66 @@ + + + + + + + + + O H C N Na + + H O N NO2 + + 2.165 + + + + + + + + + + + H:1 + + + + + 2.165 + + + + + O:1 + + + + + 2.165 + + + + + N:1 + + + + + 2.165 + + + + + O:2 N:1 + + + + + + + 2.165 + + + + + diff --git a/test_problems/statmech_test_poly/output_blessed.txt b/test_problems/statmech_test_poly/output_blessed.txt new file mode 100644 index 000000000..1fc388107 --- /dev/null +++ b/test_problems/statmech_test_poly/output_blessed.txt @@ -0,0 +1,11 @@ + + +************************************************ + Cantera Error! +************************************************ + + +Procedure: installStatMechThermoFromXML +Error: Expected no coeff: this is not a polynomial representation + + diff --git a/test_problems/statmech_test_poly/statmech_test_poly.cpp b/test_problems/statmech_test_poly/statmech_test_poly.cpp new file mode 100644 index 000000000..640613c78 --- /dev/null +++ b/test_problems/statmech_test_poly/statmech_test_poly.cpp @@ -0,0 +1,50 @@ +/** + * @file statmech + * test problem for statistical mechanics in cantera + */ + +// Example +// +// Test case for the statistical mechanics in cantera +// + +#include "cantera/transport.h" +#include "cantera/IdealGasMix.h" +#include "cantera/equil/equil.h" + +using namespace std; +using namespace Cantera; + +int main(int argc, char** argv) +{ + + try { + int k; + IdealGasMix g("test_stat_err.xml"); + int nsp = g.nSpecies(); + double pres = 1.0E5; + + vector_fp Xset(nsp, 0.0); + Xset[0] = 0.5 ; + Xset[1] = 0.5; + + g.setState_TPX(1500.0, pres, DATA_PTR(Xset)); + equilibrate(g, "TP", -1); + + vector_fp cp_R(nsp, 0.0); + g.getCp_R(DATA_PTR(cp_R)); + + } catch (CanteraError) { + // we wanted to catch an error here for incorrectly trying to use poly methods + // for the statmech species data, so exit successfully + showErrors(cout); + + // Mark it zero! + return 0; + } + + cout << "ERROR" << std::endl; + // something is wrong: we were suppose to catch an error here: paradox! + return 1; + +} diff --git a/test_problems/statmech_test_poly/test_stat_err.xml b/test_problems/statmech_test_poly/test_stat_err.xml new file mode 100644 index 000000000..5af385f26 --- /dev/null +++ b/test_problems/statmech_test_poly/test_stat_err.xml @@ -0,0 +1,47 @@ + + + + + + + + + O H C N Na Cl + + H O + + 2.165 + + + + + + + + + + + H:1 + + + + + 2.165 + + + + + O:1 + + + + 2.344331120E+00, 7.980520750E-03, -1.947815100E-05, 2.015720940E-08, + -7.376117610E-12, -9.179351730E+02, 6.830102380E-01 + + + 2.165 + + + + + diff --git a/test_problems/statmech_transport/output_blessed.txt b/test_problems/statmech_transport/output_blessed.txt new file mode 100644 index 000000000..012ea92e8 --- /dev/null +++ b/test_problems/statmech_transport/output_blessed.txt @@ -0,0 +1 @@ +here diff --git a/test_problems/statmech_transport/statmech_transport.cpp b/test_problems/statmech_transport/statmech_transport.cpp new file mode 100644 index 000000000..51ed0f69f --- /dev/null +++ b/test_problems/statmech_transport/statmech_transport.cpp @@ -0,0 +1,103 @@ +/** + * @file statmech + * test problem for statistical mechanics in cantera + */ + +// Example +// +// Test case for the statistical mechanics in cantera +// + +#include "cantera/transport.h" +#include "cantera/IdealGasMix.h" +#include "cantera/equil/equil.h" + +using namespace std; +using namespace Cantera; + +int main(int argc, char** argv) +{ + + try { + int k; + IdealGasMix g("test_stat_trans.xml", "example"); + int nsp = g.nSpecies(); + double pres = 1.0E5; + + vector_fp Xset(nsp, 0.0); + Xset[0] = 0.5 ; + Xset[1] = 0.5; + + g.setState_TPX(1500.0, pres, DATA_PTR(Xset)); + equilibrate(g, "TP", -1); + + // init pecos transport + int log_level = 0; + Transport* tran = newTransportMgr("Pecos", &g, log_level=0); + PecosTransport* tranMix = dynamic_cast(tran); + + cout << "here" << std::endl; + + + vector_fp cp_R(nsp, 0.0); + g.getCp_R(DATA_PTR(cp_R)); + + //for(int i=0;i= tol) { + double diff = cp_R[3]-sol; + std::cout << "Error for Species NO2!\n"; + std::cout << "Diff was: " << diff << "\n"; + return 1; + } + + } catch (CanteraError) { + showErrors(cout); + return 1; + } + + // Mark it zero! + return 0; + +} diff --git a/test_problems/statmech_transport/test_stat_trans.xml b/test_problems/statmech_transport/test_stat_trans.xml new file mode 100644 index 000000000..91673d957 --- /dev/null +++ b/test_problems/statmech_transport/test_stat_trans.xml @@ -0,0 +1,100 @@ + + + + + + + + + O H C N Na + + H O N NO2 + + 2.165 + + + + + + + + + + + H:1 + + + + + 2.165 + + atom + 80.000 + 2.750 + 0.000 + 0.000 + 0.000 + + + + + + O:1 + + + + + 2.165 + + atom + 80.000 + 2.750 + 0.000 + 0.000 + 0.000 + + + + + + N:1 + + + + + 2.165 + + atom + 80.000 + 2.750 + 0.000 + 0.000 + 0.000 + + + + + + O:2 N:1 + + + + + + + 2.165 + + + atom + 80.000 + 2.750 + 0.000 + 0.000 + 0.000 + + + + + + +