research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Added back Statmech and PecosTransport.

Browse source 
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.
This commit is contained in:
Harry Moffat 2014-03-15 02:57:00 +00:00
parent 3d5aa42e77
commit 39651964a7
34 changed files with 3152 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
include/cantera/thermo/GeneralSpeciesThermo.h
View file

@ -12,6 +12,7 @@
#include "SpeciesThermoMgr.h" #include "SpeciesThermoMgr.h"
#include "NasaPoly1.h" #include "NasaPoly1.h"
#include "Nasa9Poly1.h" #include "Nasa9Poly1.h"
#include "StatMech.h"
#include "speciesThermoTypes.h" #include "speciesThermoTypes.h"
namespace Cantera namespace Cantera

59
include/cantera/thermo/IdealGasPhase.h
View file

@ -439,6 +439,65 @@ public:
*/ */
virtual doublereal cv_mole() const; 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 //! @name Mechanical Equation of State
//! @{ //! @{

171
include/cantera/thermo/StatMech.h Normal file
View file

@ -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<std::string,species*> name_map;
};
}
#endif

4
include/cantera/thermo/speciesThermoTypes.h
View file

@ -58,6 +58,10 @@
//! This is implemented in the class Nasa9PolyMultiTempRegion in Nasa9Poly1MultiTempRegion //! This is implemented in the class Nasa9PolyMultiTempRegion in Nasa9Poly1MultiTempRegion
#define NASA9MULTITEMP 513 #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. //! Surface Adsorbate Model for a species on a surface.
//! This is implemented in the class Adsorbate. //! This is implemented in the class Adsorbate.
#define ADSORBATE 1024 #define ADSORBATE 1024

1
include/cantera/transport.h
View file

@ -13,5 +13,6 @@
#include "transport/DustyGasTransport.h" #include "transport/DustyGasTransport.h"
#include "transport/MultiTransport.h" #include "transport/MultiTransport.h"
#include "transport/MixTransport.h" #include "transport/MixTransport.h"
#include "transport/PecosTransport.h"
#include "transport/LiquidTransport.h" #include "transport/LiquidTransport.h"
#endif #endif

314
include/cantera/transport/PecosTransport.h Normal file
View file

@ -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<vector_fp> m_visccoeffs;
std::vector<vector_fp> m_condcoeffs;
std::vector<vector_fp> 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<std::vector<int> > m_poly;
std::vector<vector_fp> m_astar_poly;
std::vector<vector_fp> m_bstar_poly;
std::vector<vector_fp> m_cstar_poly;
std::vector<vector_fp> 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

1
include/cantera/transport/TransportBase.h
View file

@ -54,6 +54,7 @@ const int cAqueousTransport = 750;
const int cSimpleTransport = 770; const int cSimpleTransport = 770;
const int cRadiativeTransport = 800; const int cRadiativeTransport = 800;
const int cWaterTransport = 721; const int cWaterTransport = 721;
const int cPecosTransport = 900;
//! \endcond //! \endcond
// forward reference // forward reference

5
src/thermo/GeneralSpeciesThermo.cpp
View file

@ -145,6 +145,11 @@ void GeneralSpeciesThermo::install(const std::string& name,
refPressure_, c); refPressure_, c);
break; break;
case STAT:
m_sp[index] = new StatMech(index, minTemp_, maxTemp_,
refPressure_, c, name);
break;
case ADSORBATE: case ADSORBATE:
m_sp[index] = new Adsorbate(index, minTemp_, maxTemp_, m_sp[index] = new Adsorbate(index, minTemp_, maxTemp_,
refPressure_, c); refPressure_, c);

55
src/thermo/IdealGasPhase.cpp
View file

@ -94,6 +94,61 @@ doublereal IdealGasPhase::cv_mole() const
return cp_mole() - GasConstant; 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 doublereal IdealGasPhase::standardConcentration(size_t k) const
{ {
double p = pressure(); double p = pressure();

53
src/thermo/SpeciesThermoFactory.cpp
View file

@ -17,6 +17,7 @@ using namespace std;
#include "cantera/thermo/Mu0Poly.h" #include "cantera/thermo/Mu0Poly.h"
#include "Nasa9PolyMultiTempRegion.h" #include "Nasa9PolyMultiTempRegion.h"
#include "cantera/thermo/Nasa9Poly1.h" #include "cantera/thermo/Nasa9Poly1.h"
#include "cantera/thermo/StatMech.h"
#include "cantera/thermo/AdsorbateThermo.h" #include "cantera/thermo/AdsorbateThermo.h"
#include "cantera/thermo/SpeciesThermoMgr.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<XML_Node*>& tp)
{
const XML_Node* fptr = tp[0];
int nRegTmp = tp.size();
vector_fp cPoly;
std::vector<StatMech*> 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. //! 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. * 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); installMu0ThermoFromXML(speciesNode["name"], spthermo, k, f);
} else if (f->name() == "NASA9") { } 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") { } else if (f->name() == "adsorbate") {
installAdsorbateThermoFromXML(speciesNode["name"], spthermo, k, *f); installAdsorbateThermoFromXML(speciesNode["name"], spthermo, k, *f);
} else { } else {
@ -678,6 +727,8 @@ void SpeciesThermoFactory::installThermoForSpecies
installNasaThermoFromXML(speciesNode["name"], spthermo, k, f0, f1); installNasaThermoFromXML(speciesNode["name"], spthermo, k, f0, f1);
} else if (f0->name() == "Shomate" && f1->name() == "Shomate") { } else if (f0->name() == "Shomate" && f1->name() == "Shomate") {
installShomateThermoFromXML(speciesNode["name"], spthermo, k, f0, f1); 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") { } else if (f0->name() == "NASA9" && f1->name() == "NASA9") {
installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
} else { } else {
@ -688,6 +739,8 @@ void SpeciesThermoFactory::installThermoForSpecies
const XML_Node* f0 = tp[0]; const XML_Node* f0 = tp[0];
if (f0->name() == "NASA9") { if (f0->name() == "NASA9") {
installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp); installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
} else if (f0->name() == "StatMech") {
installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp);
} else { } else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"], throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
"multiple"); "multiple");

655
src/thermo/StatMech.cpp Normal file
View file

@ -0,0 +1,655 @@
/**
* @file StatMech.cpp
* \link Cantera::SpeciesThermoInterpType SpeciesThermoInterpType\endlink
*/
// Copyright 2007 Sandia National Laboratories
#include "cantera/thermo/StatMech.h"
#include <iostream>
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<std::string> 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; i<name_map[SS[ii]]->nvib; 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<std::string,species*>::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)
{
}
}

592
src/transport/PecosTransport.cpp Executable file
View file

@ -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 <sstream>
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; k<m_nsp; k++) {
double sum1 = 0.0;
double sum2 = 0.0;
for (int i=0; i<m_nsp; i++) {
if (i==k) {
continue;
}
sum1 += m_molefracs[i] / m_bdiff(k,i);
sum2 += m_molefracs[i] * m_mw[i] / m_bdiff(k,i);
}
sum1 *= p;
sum2 *= p * m_molefracs[k] / (mmw - m_mw[k]*m_molefracs[k]);
d[k] = 1.0 / (sum1 + sum2);
}
}
}
/**
* @internal This is called whenever a transport property is
* requested from ThermoSubstance if the temperature has changed
* since the last call to update_T.
*/
void PecosTransport::update_T()
{
doublereal t = m_thermo->temperature();
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;
}
}

8
src/transport/TransportFactory.cpp
View file

@ -8,6 +8,7 @@
// known transport models // known transport models
#include "cantera/transport/MultiTransport.h" #include "cantera/transport/MultiTransport.h"
#include "cantera/transport/PecosTransport.h"
#include "cantera/transport/MixTransport.h" #include "cantera/transport/MixTransport.h"
#include "cantera/transport/SolidTransport.h" #include "cantera/transport/SolidTransport.h"
#include "cantera/transport/DustyGasTransport.h" #include "cantera/transport/DustyGasTransport.h"
@ -175,6 +176,7 @@ TransportFactory::TransportFactory() :
m_models["Aqueous"] = cAqueousTransport; m_models["Aqueous"] = cAqueousTransport;
m_models["Simple"] = cSimpleTransport; m_models["Simple"] = cSimpleTransport;
m_models["User"] = cUserTransport; m_models["User"] = cUserTransport;
m_models["Pecos"] = cPecosTransport;
m_models["None"] = None; m_models["None"] = None;
//m_models["Radiative"] = cRadiative; //m_models["Radiative"] = cRadiative;
for (map<string, int>::iterator iter = m_models.begin(); for (map<string, int>::iterator iter = m_models.begin();
@ -341,7 +343,13 @@ Transport* TransportFactory::newTransport(const std::string& transportModel,
tr = new MixTransport; tr = new MixTransport;
initTransport(tr, phase, CK_Mode, log_level); initTransport(tr, phase, CK_Mode, log_level);
break; break;
// adding pecos transport model 2/13/12
case cPecosTransport:
tr = new PecosTransport;
initTransport(tr, phase, 0, log_level);
break;
case cSolidTransport: case cSolidTransport:
tr = new SolidTransport; tr = new SolidTransport;
initSolidTransport(tr, phase, log_level); initSolidTransport(tr, phase, log_level);
tr->setThermo(*phase); tr->setThermo(*phase);

252
test_problems/PecosTransport/PecosTransport.cpp Normal file
View file

@ -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 <cstdio>
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<PecosTransport*>(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;
}
/***********************************************************/

176
test_problems/PecosTransport/output_blessed.txt Normal file
View file

@ -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!.

36
test_problems/PecosTransport/runtest Executable file
View file

@ -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

11
test_problems/SConscript
View file

@ -297,6 +297,7 @@ if haveConverters:
Test('negA-cti', 'negATest', negA, 'negATest_blessed.out', Test('negA-cti', 'negATest', negA, 'negATest_blessed.out',
options='noxNeg.cti', artifacts=negA_name) options='noxNeg.cti', artifacts=negA_name)
CompileAndTest('pecosTransport', 'PecosTransport', 'pecosTransport', 'output_blessed.txt')
CompileAndTest('printUtil', 'printUtilUnitTest', 'pUtest', 'output_blessed.txt') CompileAndTest('printUtil', 'printUtilUnitTest', 'pUtest', 'output_blessed.txt')
CompileAndTest('pureFluid', 'pureFluidTest', 'testPureWater', 'output_blessed.txt') CompileAndTest('pureFluid', 'pureFluidTest', 'testPureWater', 'output_blessed.txt')
if haveConverters: if haveConverters:
@ -307,6 +308,16 @@ CompileAndTest('simpleTransport', 'simpleTransport', 'simpleTransport',
'output_blessed.txt') 'output_blessed.txt')
CompileAndTest('stoichSolidKinetics', 'stoichSolidKinetics', CompileAndTest('stoichSolidKinetics', 'stoichSolidKinetics',
'stoichSolidKinetics', 'output_blessed.txt') '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('surfkin', 'surfkin', 'surfdemo', 'output_blessed.txt')
CompileAndTest('surfSolver', 'surfSolverTest', 'surfaceSolver', None, CompileAndTest('surfSolver', 'surfSolverTest', 'surfaceSolver', None,
arguments='haca2.xml', arguments='haca2.xml',

23
test_problems/statmech/output_blessed.txt Normal file
View file

@ -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!.

15
test_problems/statmech/runtest_stat Executable file
View file

@ -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 $?

57
test_problems/statmech/statmech_test_Fe.cpp Normal file
View file

@ -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<nsp; i++) {
std::cout << cp_R[i] << std::endl;
}
// error check
if (cp_R[0] != 0) {
std::cout << "Error for monotomic Species!\n";
return 1;
}
} catch (CanteraError) {
// need to get error here because of loading Fe in input file, when
// no Fe exists in the species information table, in statmech.cpp
showErrors(cout);
return 0;
}
// Mark it zero!
return 1;
}

44
test_problems/statmech/test_stat_Fe.xml Normal file
View file

@ -0,0 +1,44 @@
<?xml version="1.0"?>
<ctml>
<validate reactions="yes" species="yes"/>
<!-- phase H -->
<phase dim="3" id="H">
<elementArray datasrc="elements.xml">
O H C N Na Cl Fe
</elementArray>
<speciesArray datasrc="#species_test"> H Fe</speciesArray>
<thermo model="IdealGas">
<density units="g/cm3">2.165</density>
</thermo>
<transport model="None"/>
<kinetics model="none"/>
</phase>
<!-- species definitions -->
<speciesData id="species_test">
<!-- species H -->
<species name="H">
<atomArray> H:1 </atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="2000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species Fe -->
<species name="Fe">
<atomArray> Fe:1 </atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="2000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
</speciesData>
</ctml>

4
test_problems/statmech_properties/output_blessed.txt Normal file
View file

@ -0,0 +1,4 @@
2.5
2.5
2.5
6.6558

61
test_problems/statmech_properties/statmech_properties.cpp Normal file
View file

@ -0,0 +1,61 @@
/**
* @file statmech
* test problem for statistical mechanics in cantera
*/
// Example
//
// Test case for the statistical mechanics in cantera
//
#include <iostream>
#include <string>
#include <vector>
#include <string>
#include <iomanip>
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;
}

66
test_problems/statmech_properties/test_stat.xml Normal file
View file

@ -0,0 +1,66 @@
<?xml version="1.0"?>
<ctml>
<validate reactions="yes" species="yes"/>
<!-- phase H -->
<phase dim="3" id="H">
<elementArray datasrc="elements.xml">
O H C N Na
</elementArray>
<speciesArray datasrc="#species_test"> H O N NO2</speciesArray>
<thermo model="IdealGas">
<density units="g/cm3">2.165</density>
</thermo>
<transport model="None"/>
<kinetics model="none"/>
</phase>
<!-- species definitions -->
<speciesData id="species_test">
<!-- species H -->
<species name="H">
<atomArray> H:1 </atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species O -->
<species name="O">
<atomArray>O:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species N -->
<species name="N">
<atomArray>N:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species NO2 -->
<species name="NO2">
<atomArray>O:2 N:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
</speciesData>
</ctml>

0
test_problems/statmech_test/output_blessed.txt Normal file
View file

15
test_problems/statmech_test/runtest_stat Executable file
View file

@ -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 $?

95
test_problems/statmech_test/statmech_test.cpp Normal file
View file

@ -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<nsp;i++)
//{
// std::cout.precision(10);
// std::cout << cp_R[i] << std::endl;
// }
// error check-- exactly 2.5 for atoms
if (cp_R[0] != 2.5) {
std::cout << "Error for monotomic Species!\n";
return 1;
}
// error check: analytical result is more complicated for
// molecules. One species should suffice, lets try NO2, with
// three vibrational modes:
/// theta[0]: 1.07900e3
/// theta[1]: 1.90000e3
/// theta[2]: 2.32700e3
// at T = 1500
//
// This is precisely: 6.655804161 (e.g. 5/2 + 2 + 3.1558..)
//
double theta[3];
theta[0] = 1.07900e3;
theta[1] = 1.90000e3;
theta[2] = 2.32700e3;
double T;
T = 1500.0;
double denom;
double ctr = 0.0;
double GasConstant = 1.0;
for (int i = 0; i < 3; i++) {
denom = exp(2*theta[i]/T) - 2* exp(theta[i]/T) + 1;
ctr += GasConstant * theta[i] * (theta[i] * exp(theta[i]/T)/(T*T))/ (denom);
//std::cout << "survey says: " << ctr << " and denom is: " << denom << std::endl;
}
//std::cout << "survey says: " << ctr << " and denom is: " << denom << std::endl;
double sol = ctr + 5/2 + 2;
double tol = 1e-9;
if (abs(cp_R[3] - sol) >= 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;
}

66
test_problems/statmech_test/test_stat.xml Normal file
View file

@ -0,0 +1,66 @@
<?xml version="1.0"?>
<ctml>
<validate reactions="yes" species="yes"/>
<!-- phase H -->
<phase dim="3" id="H">
<elementArray datasrc="elements.xml">
O H C N Na
</elementArray>
<speciesArray datasrc="#species_test"> H O N NO2</speciesArray>
<thermo model="IdealGas">
<density units="g/cm3">2.165</density>
</thermo>
<transport model="None"/>
<kinetics model="none"/>
</phase>
<!-- species definitions -->
<speciesData id="species_test">
<!-- species H -->
<species name="H">
<atomArray> H:1 </atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species O -->
<species name="O">
<atomArray>O:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species N -->
<species name="N">
<atomArray>N:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species NO2 -->
<species name="NO2">
<atomArray>O:2 N:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
</speciesData>
</ctml>

11
test_problems/statmech_test_poly/output_blessed.txt Normal file
View file

@ -0,0 +1,11 @@
************************************************
Cantera Error!
************************************************
Procedure: installStatMechThermoFromXML
Error: Expected no coeff: this is not a polynomial representation

50
test_problems/statmech_test_poly/statmech_test_poly.cpp Normal file
View file

@ -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;
}

47
test_problems/statmech_test_poly/test_stat_err.xml Normal file
View file

@ -0,0 +1,47 @@
<?xml version="1.0"?>
<ctml>
<validate reactions="yes" species="yes"/>
<!-- phase H -->
<phase dim="3" id="H">
<elementArray datasrc="elements.xml">
O H C N Na Cl
</elementArray>
<speciesArray datasrc="#species_test"> H O</speciesArray>
<thermo model="IdealGas">
<density units="g/cm3">2.165</density>
</thermo>
<transport model="None"/>
<kinetics model="none"/>
</phase>
<!-- species definitions -->
<speciesData id="species_test">
<!-- species H -->
<species name="H">
<atomArray> H:1 </atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="1000.0" Tmin="0.1">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
<!-- species O -->
<species name="O">
<atomArray> O:1 </atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="1000.0" Tmin="0.1">
<floatArray name="coeffs" size="7">
2.344331120E+00, 7.980520750E-03, -1.947815100E-05, 2.015720940E-08,
-7.376117610E-12, -9.179351730E+02, 6.830102380E-01</floatArray>
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
</species>
</speciesData>
</ctml>

1
test_problems/statmech_transport/output_blessed.txt Normal file
View file

@ -0,0 +1 @@
here

103
test_problems/statmech_transport/statmech_transport.cpp Normal file
View file

@ -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<PecosTransport*>(tran);
cout << "here" << std::endl;
vector_fp cp_R(nsp, 0.0);
g.getCp_R(DATA_PTR(cp_R));
//for(int i=0;i<nsp;i++)
//{
// std::cout.precision(10);
// std::cout << cp_R[i] << std::endl;
// }
// error check-- exactly 2.5 for atoms
if (cp_R[0] != 2.5) {
std::cout << "Error for monotomic Species!\n";
return 1;
}
// error check: analytical result is more complicated for
// molecules. One species should suffice, lets try NO2, with
// three vibrational modes:
/// theta[0]: 1.07900e3
/// theta[1]: 1.90000e3
/// theta[2]: 2.32700e3
// at T = 1500
//
// This is precisely: 6.655804161 (e.g. 5/2 + 2 + 3.1558..)
//
double theta[3];
theta[0] = 1.07900e3;
theta[1] = 1.90000e3;
theta[2] = 2.32700e3;
double T;
T = 1500.0;
double denom;
double ctr = 0.0;
double GasConstant = 1.0;
for (int i = 0; i < 3; i++) {
denom = exp(2*theta[i]/T) - 2* exp(theta[i]/T) + 1;
ctr += GasConstant * theta[i] * (theta[i] * exp(theta[i]/T)/(T*T))/ (denom);
//std::cout << "survey says: " << ctr << " and denom is: " << denom << std::endl;
}
//std::cout << "survey says: " << ctr << " and denom is: " << denom << std::endl;
double sol = ctr + 5/2 + 2;
double tol = 1e-9;
if (abs(cp_R[3] - sol) >= 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;
}

100
test_problems/statmech_transport/test_stat_trans.xml Normal file
View file

@ -0,0 +1,100 @@
<?xml version="1.0"?>
<ctml>
<validate reactions="yes" species="yes"/>
<!-- phase H -->
<phase dim="3" id="example">
<elementArray datasrc="elements.xml">
O H C N Na
</elementArray>
<speciesArray datasrc="#species_test"> H O N NO2</speciesArray>
<thermo model="IdealGas">
<density units="g/cm3">2.165</density>
</thermo>
<kinetics model="none"/>
<transport model="Pecos"/>
</phase>
<!-- species definitions -->
<speciesData id="species_test">
<!-- species H -->
<species name="H">
<atomArray> H:1 </atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
<transport model="Pecos">
<string title="geometry">atom</string>
<LJ_welldepth units="K">80.000</LJ_welldepth>
<LJ_diameter units="A">2.750</LJ_diameter>
<dipoleMoment units="Debye">0.000</dipoleMoment>
<polarizability units="A3">0.000</polarizability>
<rotRelax>0.000</rotRelax>
</transport>
</species>
<!-- species O -->
<species name="O">
<atomArray>O:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
<transport model="Pecos">
<string title="geometry">atom</string>
<LJ_welldepth units="K">80.000</LJ_welldepth>
<LJ_diameter units="A">2.750</LJ_diameter>
<dipoleMoment units="Debye">0.000</dipoleMoment>
<polarizability units="A3">0.000</polarizability>
<rotRelax>0.000</rotRelax>
</transport>
</species>
<!-- species N -->
<species name="N">
<atomArray>N:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
<transport model="Pecos">
<string title="geometry">atom</string>
<LJ_welldepth units="K">80.000</LJ_welldepth>
<LJ_diameter units="A">2.750</LJ_diameter>
<dipoleMoment units="Debye">0.000</dipoleMoment>
<polarizability units="A3">0.000</polarizability>
<rotRelax>0.000</rotRelax>
</transport>
</species>
<!-- species NO2 -->
<species name="NO2">
<atomArray>O:2 N:1</atomArray>
<thermo>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
<StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
</StatMech>
</thermo>
<density units="g/cm3">2.165</density>
<transport model="Pecos">
<string title="geometry">atom</string>
<LJ_welldepth units="K">80.000</LJ_welldepth>
<LJ_diameter units="A">2.750</LJ_diameter>
<dipoleMoment units="Debye">0.000</dipoleMoment>
<polarizability units="A3">0.000</polarizability>
<rotRelax>0.000</rotRelax>
</transport>
</species>
</speciesData>
</ctml>