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

merging into cantera 2.0 the pecos autotools trunk

Browse source
This commit is contained in:
Nicholas Malaya 2012-07-17 21:12:16 +00:00
parent bc0a6dfabd
commit d9381b55da
18 changed files with 2600 additions and 264 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

20
bin/exp3to2.sh Executable file
View file

@ -0,0 +1,20 @@
#/bin/sh
#
# This sed script replaces 3 character exponents
# starting with 0 with 2 characters
# e-0xx -> e-xx
# e0xx -> exx
# E-0xx -> E-xx
# E0xx -> Exx
# where
# x is a digit
#
# It takes one argument, the file to be operated on.
# And, it writes to standard out. It may be used to do a
# replacement in place.
#
cp $1 .exp.txt
cat .exp.txt | sed 's/\([eE]-\)\(0\)\([0-9][0-9]\)/\1\3/g' | \
sed 's/\([eE]\)\(0\)\([0-9][0-9]\)/\1\3/g' | \
sed 's/\([eE]+\)\(0\)\([0-9][0-9]\)/\1\3/g'
rm .exp.txt

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

@ -13,26 +13,27 @@
#include "SpeciesThermoMgr.h"
#include "NasaPoly1.h"
#include "Nasa9Poly1.h"
#include "StatMech.h"
#include "speciesThermoTypes.h"
namespace Cantera
{
//! A species thermodynamic property manager for a phase.
/*!
* This is a general manager that can handle a wide variety
* of species thermodynamic polynomials for individual species.
* It is slow, however, because it recomputes the functions of
* temperature needed for each species. What it does is to create
* a vector of SpeciesThermoInterpType objects.
*
* @ingroup mgrsrefcalc
*/
class GeneralSpeciesThermo : public SpeciesThermo
{
//! A species thermodynamic property manager for a phase.
/*!
* This is a general manager that can handle a wide variety
* of species thermodynamic polynomials for individual species.
* It is slow, however, because it recomputes the functions of
* temperature needed for each species. What it does is to create
* a vector of SpeciesThermoInterpType objects.
*
* @ingroup mgrsrefcalc
*/
class GeneralSpeciesThermo : public SpeciesThermo
{
public:
public:
//! Constructor
GeneralSpeciesThermo();
@ -213,7 +214,7 @@ public:
#endif
private:
private:
//! Provide the SpeciesthermoInterpType object
/*!
* provide access to the SpeciesThermoInterpType object.
@ -225,7 +226,7 @@ private:
*/
SpeciesThermoInterpType* provideSTIT(size_t k);
protected:
protected:
/**
* This is the main unknown in the object. It is
@ -260,7 +261,7 @@ protected:
friend class VPSSMgr;
};
};
}

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

@ -424,14 +424,64 @@ public:
* property manager.
* @see SpeciesThermo
*/
virtual doublereal cp_mole() const;
virtual doublereal cp_mole() const;
/**
* Molar heat capacity at constant volume. Units: J/kmol/K.
* For an ideal gas mixture,
* \f[ \hat c_v = \hat c_p - \hat R. \f]
*/
virtual doublereal cv_mole() const;
/**
* Molar heat capacity at constant volume. Units: J/kmol/K.
* For an ideal gas mixture,
* \f[ \hat c_v = \hat c_p - \hat R. \f]
*/
virtual doublereal cv_mole() const;
/**
* @returns species translational/rotational specific heat at
* constant volume. Inferred from the species gas
* constant and number of translational/rotational
* degrees of freedom. The translational/rotational
* modes are assumed to be fully populated, and are
* given by
* \f[
* C^{tr}_{v,s} \equiv \frac{\partial e^{tr}_s}{\partial T} = \frac{5}{2} R_s
* \f]
* for diatomic molecules and
* \f[
* C^{tr}_{v,s} \equiv \frac{\partial e^{tr}_s}{\partial T} = \frac{3}{2} R_s
* \f]
* for atoms.
*/
virtual doublereal cv_tr(doublereal ) const;
/**
* @returns species translational specific heat at constant volume.
* Since the translational modes are assumed to be fully populated
* this is simply
* \f[
* C^{trans}_{v,s} \equiv \frac{\partial e^{trans}_s}{\partial T} = \frac{3}{2} R_s
* \f]
*/
virtual doublereal cv_trans() const;
/**
* @returns species rotational specific heat at constant volume.
* By convention, we lump the translational/rotational components
* \f[
* C^{tr}_{v,s} \equiv C^{trans}_{v,s} + C^{rot}_{v,s}
* \f]
* so then
* \f[
* C^{rot}_{v,s} \equiv C^{tr}_{v,s} - C^{trans}_{v,s}
* \f]
*/
virtual doublereal cv_rot(double atomicity) const;
/**
* @returns species vibrational specific heat at
* constant volume. This is defined as
* \f[
* C^{vib}_{v,s} \equiv \frac{\partial e^{vib}_s}{\partial T_V} = \frac{R_s \theta_{vs}^2 \exp\left(\theta_{vs}/T_V\right)}{\left[\left(\exp\left(\theta_{vs}/T_V\right)-1\right)T_V\right]^2}
* \f]
*/
virtual doublereal cv_vib(int k, doublereal T) const;
//@}

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

@ -0,0 +1,220 @@
/**
* @file StatMech.h
* Header for a single-species standard state object derived
* from
*/
/*
* Copywrite (2006) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#ifndef CT_STATMECH_H
#define CT_STATMECH_H
/*
* $Revision: 279 $
* $Date: 2009-12-05 13:08:43 -0600 (Sat, 05 Dec 2009) $
*/
#include "cantera/base/global.h"
#include "SpeciesThermoInterpType.h"
#include "SpeciesThermoMgr.h"
#include <string>
#include <map>
namespace Cantera {
//!
/*!
* @ingroup spthermo
*/
class StatMech : public SpeciesThermoInterpType {
public:
//! Empty constructor
StatMech();
//! constructor used in templated instantiations
/*!
* @param n Species index
* @param tlow Minimum temperature
* @param thigh Maximum temperature
* @param pref reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
StatMech(int n, doublereal tlow, doublereal thigh, doublereal pref,
const doublereal* coeffs, std::string my_name);
//! copy constructor
/*!
* @param b object to be copied
*/
StatMech(const StatMech& b);
//! assignment operator
/*!
* @param b object to be copied
*/
StatMech& operator=(const StatMech& b);
//! Destructor
virtual ~StatMech();
//! duplicator
virtual SpeciesThermoInterpType *
duplMyselfAsSpeciesThermoInterpType() const;
//! Returns the minimum temperature that the thermo
//! parameterization is valid
virtual doublereal minTemp() const;
//! Returns the maximum temperature that the thermo
//! parameterization is valid
virtual doublereal maxTemp() const;
//! Returns the reference pressure (Pa)
virtual doublereal refPressure() const;
//! Returns an integer representing the type of parameterization
virtual int reportType() const;
//! Returns an integer representing the species index
virtual size_t speciesIndex() const;
//! Build a series of maps for the properties needed for species
int buildmap();
//! Update the properties for this species, given a temperature polynomial
/*!
* This method is called with a pointer to an array containing the functions of
* temperature needed by this parameterization, and three pointers to arrays where the
* computed property values should be written. This method updates only one value in
* each array.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param tt vector of temperature polynomials
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
virtual void updateProperties(const doublereal* tt,
doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const;
//! Compute the reference-state property of one species
/*!
* Given temperature T in K, this method updates the values of
* the non-dimensional heat capacity at constant pressure,
* enthalpy, and entropy, at the reference pressure, Pref
* of one of the species. The species index is used
* to reference into the cp_R, h_RT, and s_R arrays.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param temp Temperature (Kelvin)
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
virtual void updatePropertiesTemp(const doublereal temp,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const;
//!This utility function reports back the type of
//! parameterization and all of the parameters for the
//! species, index.
/*!
* All parameters are output variables
*
* @param n Species index
* @param type Integer type of the standard type
* @param tlow output - Minimum temperature
* @param thigh output - Maximum temperature
* @param pref output - reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state. There are
* 12 of them, designed to be compatible
* with the multiple temperature formulation.
* coeffs[0] is equal to one.
* coeffs[1] is min temperature
* coeffs[2] is max temperature
* coeffs[3+i] from i =0,9 are the coefficients themselves
*/
virtual void reportParameters(size_t& n, int &type,
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const;
//! Modify parameters for the standard state
/*!
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
virtual void modifyParameters(doublereal* coeffs);
protected:
//! lowest valid temperature
doublereal m_lowT;
//! highest valid temperature
doublereal m_highT;
//! standard-state pressure
doublereal m_Pref;
//! species index
int m_index;
//! array of polynomial coefficients
vector_fp m_coeff;
std::string sp_name;
//*generic species struct that contains everything we need here
// achtung: add doxygen markup here
// achtung: convert doubles to realdoubles
struct species
{
//Nominal T-R Degrees of freedom (cv = cfs*k*T)
doublereal cfs;
// Mol. Wt. Molecular weight (kg/kmol)
doublereal mol_weight;
// number of vibrational temperatures necessary
int nvib;
// Theta_v Characteristic vibrational temperature(s) (K)
doublereal theta[5];
};
std::map<std::string,species*> name_map;
};
}
#endif

11
include/cantera/thermo/ThermoPhase.h
View file

@ -298,7 +298,16 @@ public:
return err("cv_mole");
}
/**
* @returns species vibrational specific heat at
* constant volume.
*
*/
/// Molar heat capacity at constant volume. Units: J/kmol/K.
virtual doublereal cv_vib(int, double) const {
return err("cv_vib");
}
/**
* @}
* @name Mechanical Properties

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

@ -59,6 +59,10 @@
//! This is implemented in the class Nasa9PolyMultiTempRegion in Nasa9Poly1MultiTempRegion
#define NASA9MULTITEMP 513
//! Properties derived from theoretical considerations
//! This is implemented in the class statmech in StatMech.h
#define STAT 111
//! Surface Adsorbate Model for a species on a surface.
//! This is implemented in the class Adsorbate.
#define ADSORBATE 1024

1
include/cantera/transport.h
View file

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

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

@ -0,0 +1,295 @@
/**
* @file PecosTransport.h
* Header file defining class PecosTransport
*/
/* $Author$
* $Revision$
* $Date$
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_PECOSTRAN_H
#define CT_PECOSTRAN_H
// turn off warnings under Windows
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif
// STL includes
#include <vector>
#include <string>
#include <map>
#include <numeric>
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
using namespace std;
// Cantera includes
#include "TransportBase.h"
#include "cantera/numerics/DenseMatrix.h"
namespace Cantera {
class GasTransportParams;
/**
*
* Class PecosTransport implements mixture-averaged transport
* properties for ideal gas mixtures.
*
*/
class PecosTransport : public Transport {
public:
virtual ~PecosTransport() {}
virtual int model() const { return cPecosTransport; }
//! Viscosity of the mixture
/*!
*
*/
virtual doublereal viscosity();
virtual void getSpeciesViscosities(doublereal* visc)
{ update_T(); updateViscosity_T(); copy(m_visc.begin(), m_visc.end(), visc); }
//! Return the thermal diffusion coefficients
/*!
* For this approximation, these are all zero.
*/
virtual void getThermalDiffCoeffs(doublereal* const dt);
/*! returns the mixture thermal conductivity
*
* This is computed using the lumped model,
* \f[
* k = k^{tr} + k^{ve}
* \f]
* where,
* \f[
* k^{tr}= 5/2 \mu_s C_{v,s}^{trans} + \mu_s C_{v,s}^{rot}
* \f]
* and,
* \f[
* k^{ve}= \mu_s C_{v,s}^{vib} + \mu_s C_{v,s}^{elec}
* \f]
*
*/
virtual doublereal thermalConductivity();
virtual void getBinaryDiffCoeffs(const int ld, doublereal* const d);
//! Mixture-averaged diffusion coefficients [m^2/s].
/*!
* For the single species case or the pure fluid case
* the routine returns the self-diffusion coefficient.
* This is need to avoid a Nan result in the formula
* below.
*/
virtual void getMixDiffCoeffs(doublereal* const d);
//! Returns the mixture-averaged diffusion coefficients [m^2/s].
//! These are the coefficients for calculating the molar diffusive fluxes
//! from the species mole fraction gradients, computed according to
//! Eq. 12.176 in "Chemically Reacting Flow":
//!
//! \f[ D_{km}^* = \frac{1-X_k}{\sum_{j \ne k}^K X_j/\mathcal{D}_{kj}} \f]
//!
//! @param[out] d vector of mixture-averaged diffusion coefficients for
//! each species, length m_nsp.
void getMixDiffCoeffsMole(doublereal* const d);
//! Returns the mixture-averaged diffusion coefficients [m^2/s].
//! These are the coefficients for calculating the diffusive mass fluxes
//! from the species mass fraction gradients, computed according to
//! Eq. 12.178 in "Chemically Reacting Flow":
//!
//! \f[ \frac{1}{D_{km}} = \sum_{j \ne k}^K \frac{X_j}{\mathcal{D}_{kj}} +
//! \frac{X_k}{1-Y_k} \sum_{j \ne k}^K \frac{Y_j}{\mathcal{D}_{kj}} \f]
//!
//! @param[out] d vector of mixture-averaged diffusion coefficients for
//! each species, length m_nsp.
void getMixDiffCoeffsMass(doublereal* const d);
virtual void getMobilities(doublereal* const mobil);
virtual void update_T();
virtual void update_C();
//! Get the species diffusive mass fluxes wrt to
//! the mass averaged velocity,
//! given the gradients in mole fraction and temperature
/*!
* Units for the returned fluxes are kg m-2 s-1.
*
* @param ndim Number of dimensions in the flux expressions
* @param grad_T Gradient of the temperature
* (length = ndim)
* @param ldx Leading dimension of the grad_X array
* (usually equal to m_nsp but not always)
* @param grad_X Gradients of the mole fraction
* Flat vector with the m_nsp in the inner loop.
* length = ldx * ndim
* @param ldf Leading dimension of the fluxes array
* (usually equal to m_nsp but not always)
* @param fluxes Output of the diffusive mass fluxes
* Flat vector with the m_nsp in the inner loop.
* length = ldx * ndim
*/
virtual void getSpeciesFluxes(int ndim,
const doublereal* grad_T,
int ldx,
const doublereal* grad_X,
int ldf, doublereal* fluxes);
//! Initialize the transport object
/*!
*
* Here we change all of the internal dimensions to be sufficient.
* We get the object ready to do property evaluations.
*
* @param tr Transport parameters for all of the species
* in the phase.
*
*/
virtual bool initGas( GasTransportParams& tr );
/**
*
* Reads the transport table specified (currently defaults to internal file)
*
* Reads the user-specified transport table, appending new species
* data and/or replacing default species data.
*
*/
void read_blottner_transport_table ();
friend class TransportFactory;
/**
* Return a structure containing all of the pertinent parameters
* about a species that was used to construct the Transport
* properties in this object.
*
* @param k Species number to obtain the properties from.
*/
struct GasTransportData getGasTransportData(int);
protected:
/// default constructor
PecosTransport();
private:
//! Calculate the pressure from the ideal gas law
doublereal pressure_ig() const {
return (m_thermo->molarDensity() * GasConstant *
m_thermo->temperature());
}
// mixture attributes
int m_nsp;
doublereal m_tmin, m_tmax;
vector_fp m_mw;
// polynomial fits
vector<vector_fp> m_visccoeffs;
vector<vector_fp> m_condcoeffs;
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;
vector<vector<int> > m_poly;
vector<vector_fp > m_astar_poly;
vector<vector_fp > m_bstar_poly;
vector<vector_fp > m_cstar_poly;
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();
void updateViscosity_T();
void updateCond_T();
void updateSpeciesViscosities();
void updateDiff_T();
void correctBinDiffCoeffs();
bool m_viscmix_ok;
bool m_viscwt_ok;
bool m_spvisc_ok;
bool m_diffmix_ok;
bool m_bindiff_ok;
bool m_abc_ok;
bool m_spcond_ok;
bool m_condmix_ok;
int m_mode;
DenseMatrix m_epsilon;
DenseMatrix m_diam;
DenseMatrix incl;
bool m_debug;
// specific heats
vector_fp cv_rot;
vector_fp cp_R;
vector_fp cv_int;
};
}
#endif

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

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

433
src/thermo/GeneralSpeciesThermo.cpp
View file

@ -22,106 +22,111 @@ namespace Cantera
{
/*
* Constructors
*/
GeneralSpeciesThermo::GeneralSpeciesThermo() :
/*
* Constructors
*/
GeneralSpeciesThermo::GeneralSpeciesThermo() :
SpeciesThermo(),
m_tlow_max(0.0),
m_thigh_min(1.0E30),
m_p0(OneAtm),
m_kk(0)
{
{
m_tlow_max = 0.0;
m_thigh_min = 1.0E30;
}
}
GeneralSpeciesThermo::
GeneralSpeciesThermo(const GeneralSpeciesThermo& b) :
GeneralSpeciesThermo::
GeneralSpeciesThermo(const GeneralSpeciesThermo& b) :
m_tlow_max(b.m_tlow_max),
m_thigh_min(b.m_thigh_min),
m_kk(b.m_kk)
{
{
m_sp.resize(m_kk, 0);
for (size_t k = 0; k < m_kk; k++) {
SpeciesThermoInterpType* bk = b.m_sp[k];
if (bk) {
m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType();
}
SpeciesThermoInterpType* bk = b.m_sp[k];
if (bk) {
m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType();
}
}
}
}
GeneralSpeciesThermo&
GeneralSpeciesThermo::operator=(const GeneralSpeciesThermo& b)
{
GeneralSpeciesThermo&
GeneralSpeciesThermo::operator=(const GeneralSpeciesThermo& b)
{
if (&b != this) {
m_tlow_max = b.m_tlow_max;
m_thigh_min = b.m_thigh_min;
m_tlow_max = b.m_tlow_max;
m_thigh_min = b.m_thigh_min;
for (size_t k = 0; k < m_kk; k++) {
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
delete sp;
m_sp[k] = 0;
}
}
m_kk = b.m_kk;
m_sp.resize(m_kk, 0);
for (size_t k = 0; k < m_kk; k++) {
SpeciesThermoInterpType* bk = b.m_sp[k];
if (bk) {
m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType();
}
}
for (size_t k = 0; k < m_kk; k++) {
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
delete sp;
m_sp[k] = 0;
}
}
m_kk = b.m_kk;
m_sp.resize(m_kk, 0);
for (size_t k = 0; k < m_kk; k++) {
SpeciesThermoInterpType* bk = b.m_sp[k];
if (bk) {
m_sp[k] = bk->duplMyselfAsSpeciesThermoInterpType();
}
}
}
return *this;
}
}
GeneralSpeciesThermo::~GeneralSpeciesThermo()
{
GeneralSpeciesThermo::~GeneralSpeciesThermo()
{
for (size_t k = 0; k < m_kk; k++) {
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
delete sp;
m_sp[k] = 0;
}
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
delete sp;
m_sp[k] = 0;
}
}
}
}
SpeciesThermo*
GeneralSpeciesThermo::duplMyselfAsSpeciesThermo() const
{
SpeciesThermo*
GeneralSpeciesThermo::duplMyselfAsSpeciesThermo() const
{
GeneralSpeciesThermo* gsth = new GeneralSpeciesThermo(*this);
return (SpeciesThermo*) gsth;
}
}
/*
* Install parameterization for a species.
* @param index Species index
* @param type parameterization type
* @param c coefficients. The meaning of these depends on
* the parameterization.
*/
void GeneralSpeciesThermo::install(std::string name,
size_t index,
int type,
const doublereal* c,
doublereal minTemp,
doublereal maxTemp,
doublereal refPressure)
{
/*
* Install parameterization for a species.
* @param index Species index
* @param type parameterization type
* @param c coefficients. The meaning of these depends on
* the parameterization.
*/
void GeneralSpeciesThermo::install(std::string name,
size_t index,
int type,
const doublereal* c,
doublereal minTemp,
doublereal maxTemp,
doublereal refPressure)
{
/*
* Resize the arrays if necessary, filling the empty
* slots with the zero pointer.
*/
if(minTemp <= 0.0)
{
throw CanteraError("Error in GeneralSpeciesThermo.cpp",
" Cannot take 0 tmin as input. \n\n");
}
if (index >= m_kk) {
m_sp.resize(index+1, 0);
m_kk = index+1;
m_sp.resize(index+1, 0);
m_kk = index+1;
}
//AssertThrow(m_sp[index] == 0,
// "Index position isn't null, duplication of assignment: " + int2str(index));
//int nfreq = 3;
/*
@ -130,68 +135,74 @@ void GeneralSpeciesThermo::install(std::string name,
switch (type) {
case NASA1:
m_sp[index] = new NasaPoly1(index, minTemp, maxTemp,
refPressure, c);
break;
m_sp[index] = new NasaPoly1(index, minTemp, maxTemp,
refPressure, c);
break;
case SHOMATE1:
m_sp[index] = new ShomatePoly(index, minTemp, maxTemp,
refPressure, c);
break;
m_sp[index] = new ShomatePoly(index, minTemp, maxTemp,
refPressure, c);
break;
case CONSTANT_CP:
case SIMPLE:
m_sp[index] = new ConstCpPoly(index, minTemp, maxTemp,
refPressure, c);
break;
m_sp[index] = new ConstCpPoly(index, minTemp, maxTemp,
refPressure, c);
break;
case MU0_INTERP:
m_sp[index] = new Mu0Poly(index, minTemp, maxTemp,
refPressure, c);
break;
m_sp[index] = new Mu0Poly(index, minTemp, maxTemp,
refPressure, c);
break;
case SHOMATE2:
m_sp[index] = new ShomatePoly2(index, minTemp, maxTemp,
refPressure, c);
break;
m_sp[index] = new ShomatePoly2(index, minTemp, maxTemp,
refPressure, c);
break;
case NASA2:
m_sp[index] = new NasaPoly2(index, minTemp, maxTemp,
refPressure, c);
break;
m_sp[index] = new NasaPoly2(index, minTemp, maxTemp,
refPressure, c);
break;
case STAT:
m_sp[index] = new StatMech(index, minTemp, maxTemp,
refPressure, c, name);
break;
case ADSORBATE:
m_sp[index] = new Adsorbate(index, minTemp, maxTemp,
refPressure, c);
break;
m_sp[index] = new Adsorbate(index, minTemp, maxTemp,
refPressure, c);
break;
default:
throw UnknownSpeciesThermoModel(
"GeneralSpeciesThermo::install",
"unknown species type", int2str(type));
break;
throw UnknownSpeciesThermoModel(
"GeneralSpeciesThermo::install",
"unknown species type", int2str(type));
break;
}
if (!m_sp[index]) {
cout << "Null m_sp... index = " << index << endl;
cout << "type = " << type << endl;
cout << "Null m_sp... index = " << index << endl;
cout << "type = " << type << endl;
}
m_tlow_max = max(minTemp, m_tlow_max);
m_thigh_min = min(maxTemp, m_thigh_min);
}
}
// Install a new species thermodynamic property
// parameterization for one species.
/*
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
void GeneralSpeciesThermo::install_STIT(SpeciesThermoInterpType* stit_ptr)
{
// Install a new species thermodynamic property
// parameterization for one species.
/*
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
void GeneralSpeciesThermo::install_STIT(SpeciesThermoInterpType* stit_ptr)
{
/*
* Resize the arrays if necessary, filling the empty
* slots with the zero pointer.
*/
if (!stit_ptr) {
throw CanteraError("GeneralSpeciesThermo::install_STIT",
"zero pointer");
throw CanteraError("GeneralSpeciesThermo::install_STIT",
"zero pointer");
}
size_t index = stit_ptr->speciesIndex();
if (index >= m_kk) {
m_sp.resize(index+1, 0);
m_kk = index+1;
m_sp.resize(index+1, 0);
m_kk = index+1;
}
AssertThrow(m_sp[index] == 0,
"Index position isn't null, duplication of assignment: " + int2str(index));
@ -208,161 +219,177 @@ void GeneralSpeciesThermo::install_STIT(SpeciesThermoInterpType* stit_ptr)
m_tlow_max = max(minTemp, m_tlow_max);
m_thigh_min = min(maxTemp, m_thigh_min);
}
}
void GeneralSpeciesThermo::installPDSShandler(size_t k, PDSS* PDSS_ptr,
VPSSMgr* vpssmgr_ptr)
{
void GeneralSpeciesThermo::installPDSShandler(size_t k, PDSS* PDSS_ptr,
VPSSMgr* vpssmgr_ptr)
{
STITbyPDSS* stit_ptr = new STITbyPDSS(k, vpssmgr_ptr, PDSS_ptr);
install_STIT(stit_ptr);
}
}
/**
* Update the properties for one species.
*/
void GeneralSpeciesThermo::
update_one(size_t k, doublereal t, doublereal* cp_R,
doublereal* h_RT, doublereal* s_R) const
{
/**
* Update the properties for one species.
*/
void GeneralSpeciesThermo::
update_one(size_t k, doublereal t, doublereal* cp_R,
doublereal* h_RT, doublereal* s_R) const
{
SpeciesThermoInterpType* sp_ptr = m_sp[k];
if (sp_ptr) {
sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
}
}
}
/**
* Update the properties for all species.
*/
void GeneralSpeciesThermo::
update(doublereal t, doublereal* cp_R,
doublereal* h_RT, doublereal* s_R) const
{
/**
* Update the properties for all species.
*/
void GeneralSpeciesThermo::
update(doublereal t, doublereal* cp_R,
doublereal* h_RT, doublereal* s_R) const
{
vector<SpeciesThermoInterpType*>::const_iterator _begin, _end;
_begin = m_sp.begin();
_end = m_sp.end();
SpeciesThermoInterpType* sp_ptr = 0;
for (; _begin != _end; ++_begin) {
sp_ptr = *(_begin);
if (sp_ptr) {
sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
}
// else {
// writelog("General::update: sp_ptr is NULL!\n");
//}
sp_ptr = *(_begin);
if (sp_ptr) {
sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
}
// else {
// writelog("General::update: sp_ptr is NULL!\n");
//}
}
}
}
/**
* This utility function reports the type of parameterization
* used for the species, index.
*/
int GeneralSpeciesThermo::reportType(size_t index) const
{
/**
* This utility function reports the type of parameterization
* used for the species, index.
*/
int GeneralSpeciesThermo::reportType(size_t index) const
{
SpeciesThermoInterpType* sp = m_sp[index];
if (sp) {
return sp->reportType();
return sp->reportType();
}
return -1;
}
}
/**
* This utility function reports back the type of
* parameterization and all of the parameters for the
* species, index.
* For the NASA object, there are 15 coefficients.
*/
void GeneralSpeciesThermo::
reportParams(size_t index, int& type, doublereal* const c,
doublereal& minTemp, doublereal& maxTemp, doublereal& refPressure) const
{
/**
* This utility function reports back the type of
* parameterization and all of the parameters for the
* species, index.
* For the NASA object, there are 15 coefficients.
*/
void GeneralSpeciesThermo::
reportParams(size_t index, int& type, doublereal* const c,
doublereal& minTemp, doublereal& maxTemp, doublereal& refPressure) const
{
SpeciesThermoInterpType* sp = m_sp[index];
size_t n;
if (sp) {
sp->reportParameters(n, type, minTemp, maxTemp,
refPressure, c);
if (n != index) {
throw CanteraError("GeneralSpeciesThermo::reportParams",
"Internal error encountered");
}
sp->reportParameters(n, type, minTemp, maxTemp,
refPressure, c);
if (n != index) {
throw CanteraError("GeneralSpeciesThermo::reportParams",
"Internal error encountered");
}
} else {
type = -1;
type = -1;
}
}
}
/**
* Return the lowest temperature at which the thermodynamic
* parameterization is valid. If no argument is supplied, the
* value is the one for which all species parameterizations
* are valid. Otherwise, if an integer argument is given, the
* value applies only to the species with that index.
*/
doublereal GeneralSpeciesThermo::minTemp(size_t k) const
{
// //! Modify parameters for the standard state
// /*!
// * @param index Species index
// * @param c Vector of coefficients used to set the
// * parameters for the standard state.
// */
// void GeneralSpeciesThermo::
// modifyParams(size_t index, doublereal* c)
// {
// SpeciesThermoInterpType* sp = m_sp[index];
// if (sp) {
// sp->modifyParameters(c);
// }
// }
/**
* Return the lowest temperature at which the thermodynamic
* parameterization is valid. If no argument is supplied, the
* value is the one for which all species parameterizations
* are valid. Otherwise, if an integer argument is given, the
* value applies only to the species with that index.
*/
doublereal GeneralSpeciesThermo::minTemp(size_t k) const
{
if (k == npos) {
return m_tlow_max;
return m_tlow_max;
} else {
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
return sp->minTemp();
}
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
return sp->minTemp();
}
}
return m_tlow_max;
}
}
doublereal GeneralSpeciesThermo::maxTemp(size_t k) const
{
doublereal GeneralSpeciesThermo::maxTemp(size_t k) const
{
if (k == npos) {
return m_thigh_min;
return m_thigh_min;
} else {
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
return sp->maxTemp();
}
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
return sp->maxTemp();
}
}
return m_thigh_min;
}
}
doublereal GeneralSpeciesThermo::refPressure(size_t k) const
{
doublereal GeneralSpeciesThermo::refPressure(size_t k) const
{
if (k == npos) {
return m_p0;
return m_p0;
} else {
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
return sp->refPressure();
}
SpeciesThermoInterpType* sp = m_sp[k];
if (sp) {
return sp->refPressure();
}
}
return m_p0;
}
}
SpeciesThermoInterpType* GeneralSpeciesThermo::provideSTIT(size_t k)
{
SpeciesThermoInterpType* GeneralSpeciesThermo::provideSTIT(size_t k)
{
return (m_sp[k]);
}
}
#ifdef H298MODIFY_CAPABILITY
doublereal GeneralSpeciesThermo::reportOneHf298(int k) const
{
doublereal GeneralSpeciesThermo::reportOneHf298(int k) const
{
SpeciesThermoInterpType* sp_ptr = m_sp[k];
doublereal h = -1.0;
if (sp_ptr) {
h = sp_ptr->reportHf298(0);
h = sp_ptr->reportHf298(0);
}
return h;
}
}
void GeneralSpeciesThermo::modifyOneHf298(const int k, const doublereal Hf298New)
{
void GeneralSpeciesThermo::modifyOneHf298(const int k, const doublereal Hf298New)
{
SpeciesThermoInterpType* sp_ptr = m_sp[k];
if (sp_ptr) {
sp_ptr->modifyOneHf298(k, Hf298New);
sp_ptr->modifyOneHf298(k, Hf298New);
}
}
}
#endif

92
src/thermo/IdealGasPhase.cpp
View file

@ -155,6 +155,98 @@ doublereal IdealGasPhase::cv_mole() const
return cp_mole() - GasConstant;
}
/**
* @returns species translational/rotational specific heat at
* constant volume.
*
* Either: $5/2 R_s$ or $3/2 R_s$ for molecules/atoms.
*
*/
doublereal IdealGasPhase::cv_tr (doublereal atomicity) const
{
// k is the species number
int dum = 0;
int type = 0;
doublereal c[12];
doublereal minTemp;
doublereal maxTemp;
doublereal refPressure;
m_spthermo->reportParams(dum,type,c,minTemp,maxTemp,refPressure);
if(type != 111)
{
throw CanteraError("Error in IdealGasPhase.cpp",
"cv_tr only supported for StatMech!. \n\n");
}
// see reportParameters for specific details
return c[3];
}
/**
* @returns species translational specific heat at constant volume.
*/
doublereal IdealGasPhase::cv_trans () const
{ return 1.5*GasConstant; }
/**
* @returns species rotational specific heat at constant volume.
*
*/
doublereal IdealGasPhase::cv_rot (double atom) const
{ return std::max(cv_tr(atom) - cv_trans(), 0.); }
/**
* @returns species vibrational specific heat at
* constant volume.
*
* C^{vib}_{v,s} = \frac{\partial e^{vib}_{v,s} }{\partial T}
*
* The species vibration energy ($e^{vib}_{v,s}$) is:
*
* 0: atom
*
* Diatomic:
* \f[
* \frac{R_s \theta_{v,s}}{e^{\theta_{v,s}/T}-1}
* \f]
*
* General Molecules:
* \f[
* \sum_i \frac{R_s \theta_{v,s,i}}{e^{\theta_{v,s,i}/T}-1}
* \f]
*
*/
doublereal IdealGasPhase::cv_vib (const int k, const doublereal T) const
{
// k is the species number
int dum = 0;
int type = 0;
doublereal c[12];
doublereal minTemp;
doublereal maxTemp;
doublereal refPressure;
c[0] = temperature();
m_spthermo->reportParams(dum,type,c,minTemp,maxTemp,refPressure);
// basic sanity check
if(type != 111)
{
throw CanteraError("Error in IdealGasPhase.cpp",
"cv_vib only supported for StatMech!. \n\n");
}
// see reportParameters for specific details
return c[4];
}
// Mechanical Equation of State ----------------------------
// Chemical Potentials and Activities ----------------------

27
src/thermo/Makefile.am
View file

@ -20,33 +20,10 @@ cc_sources = ConstCpPoly.cpp ConstDensityThermo.cpp DebyeHuckel.cpp \
VPSSMgrFactory.cpp VPSSMgr_ConstVol.cpp VPSSMgr_General.cpp \
VPSSMgr_IdealGas.cpp VPSSMgr_Water_ConstVol.cpp \
VPSSMgr_Water_HKFT.cpp VPStandardStateTP.cpp WaterProps.cpp \
WaterPropsIAPWS.cpp WaterPropsIAPWSphi.cpp WaterSSTP.cpp
WaterPropsIAPWS.cpp WaterPropsIAPWSphi.cpp WaterSSTP.cpp \
StatMech.cpp
#Elements.cpp Phase.cpp RedlichKisterVPSSTP.cpp \
ThermoPhase.cpp IdealGasPhase.cpp ConstDensityThermo.cpp \
SpeciesThermoFactory.cpp ConstCpPoly.cpp Nasa9Poly1.cpp \
Nasa9PolyMultiTempRegion.cpp PDSS_Water.cpp PDSS_HKFT.cpp \
Mu0Poly.cpp GeneralSpeciesThermo.cpp SurfPhase.cpp \
ThermoFactory.cpp SpeciesThermoInterpType.cpp \
VPSSMgr.cpp VPSSMgrFactory.cpp VPSSMgr_General.cpp \
IdealSolnGasVPSS.cpp MolalityVPSSTP.cpp VPStandardStateTP.cpp \
VPSSMgr_IdealGas.cpp VPSSMgr_ConstVol.cpp PDSS_ConstVol.cpp \
PDSS_IdealGas.cpp PDSS_SSVol.cpp DebyeHuckel.cpp PDSS.cpp \
WaterProps.cpp WaterPropsIAPWS.cpp WaterPropsIAPWSphi.cpp \
VPSSMgr_Water_HKFT.cpp VPSSMgr_Water_ConstVol.cpp \
PDSS_IonsFromNeutral.cpp IonsFromNeutralVPSSTP.cpp \
GibbsExcessVPSSTP.cpp LatticePhase.cpp IdealMolalSoln.cpp \
HMWSoln.cpp HMWSoln_input.cpp WaterSSTP.cpp \
MetalSHEelectrons.cpp \
IdealSolidSolnPhase.cpp LatticeSolidPhase.cpp \
SingleSpeciesTP.cpp MineralEQ3.cpp \
PseudoBinaryVPSSTP.cpp MargulesVPSSTP.cpp \
StoichSubstanceSSTP.cpp PureFluidPhase.cpp \
StoichSubstance.cpp
#PecosGasPhase.cpp
AM_CPPFLAGS = -I$(top_builddir)/include
AM_CXXFLAGS = $(AM_CPPFLAGS)

73
src/thermo/SpeciesThermoFactory.cpp
View file

@ -17,6 +17,7 @@ using namespace std;
#include "cantera/thermo/Mu0Poly.h"
#include "Nasa9PolyMultiTempRegion.h"
#include "cantera/thermo/Nasa9Poly1.h"
#include "cantera/thermo/StatMech.h"
#include "cantera/thermo/AdsorbateThermo.h"
#include "cantera/thermo/SpeciesThermoMgr.h"
@ -644,6 +645,53 @@ static void installNasa9ThermoFromXML(std::string speciesName,
}
}
/**
* Install a stat mech based property solver
* for species k into a SpeciesThermo instance.
*/
static void installStatMechThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const std::vector<XML_Node*>& tp)
{
const XML_Node * fptr = tp[0];
int nRegTmp = tp.size();
int nRegions = 0;
vector_fp cPoly;
StatMech *np_ptr = 0;
std::vector<StatMech *> regionPtrs;
doublereal tmin, tmax, pref = OneAtm;
// Loop over all of the possible temperature regions
for (int i = 0; i < nRegTmp; i++) {
fptr = tp[i];
if (fptr) {
if (fptr->name() == "StatMech") {
if (fptr->hasChild("floatArray")) {
tmin = fpValue((*fptr)["Tmin"]);
tmax = fpValue((*fptr)["Tmax"]);
if ((*fptr).hasAttrib("P0")) {
pref = fpValue((*fptr)["P0"]);
}
if ((*fptr).hasAttrib("Pref")) {
pref = fpValue((*fptr)["Pref"]);
}
getFloatArray(fptr->child("floatArray"), cPoly, false);
if (cPoly.size() != 0) {
throw CanteraError("installStatMechThermoFromXML",
"Expected no coeff: this is not a polynomial representation");
}
}
}
}
}
// set properties
tmin = 0.1;
vector_fp coeffs(1);
coeffs[0] = 0.0;
(&sp)->install(speciesName, k, STAT, &coeffs[0], tmin, tmax, pref);
}
//! Install a Adsorbate polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance.
/*!
@ -746,8 +794,11 @@ void SpeciesThermoFactory::installThermoForSpecies
} else if (f->name() == "Mu0") {
installMu0ThermoFromXML(speciesNode["name"], spthermo, k, f);
} else if (f->name() == "NASA9") {
installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
}
else if (f->name() == "StatMech") {
installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp);
}
else if (f->name() == "adsorbate") {
installAdsorbateThermoFromXML(speciesNode["name"], spthermo, k, *f);
}
@ -762,7 +813,11 @@ void SpeciesThermoFactory::installThermoForSpecies
installNasaThermoFromXML(speciesNode["name"], spthermo, k, f0, f1);
} else if (f0->name() == "Shomate" && f1->name() == "Shomate") {
installShomateThermoFromXML(speciesNode["name"], spthermo, k, f0, f1);
} else if (f0->name() == "NASA9" && f1->name() == "NASA9") {
}
else if (f0->name() == "StatMech") {
installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp);
}
else if (f0->name() == "NASA9" && f1->name() == "NASA9") {
installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
} else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
@ -772,13 +827,17 @@ void SpeciesThermoFactory::installThermoForSpecies
const XML_Node* f0 = tp[0];
if (f0->name() == "NASA9") {
installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
} else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
"multiple");
}
else if (f0->name() == "StatMech") {
installStatMechThermoFromXML(speciesNode["name"], spthermo, k, tp);
}
else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
"multiple");
}
} else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
"multiple");
throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
"multiple");
}
}
}

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

@ -0,0 +1,811 @@
/**
* @file StatMech.cpp
* \link Cantera::SpeciesThermoInterpType SpeciesThermoInterpType\endlink
*/
/* $Author: hkmoffa $
* $Revision: 279 $
* $Date: 2009-12-05 13:08:43 -0600 (Sat, 05 Dec 2009) $
*/
// Copyright 2007 Sandia National Laboratories
#include "cantera/thermo/StatMech.h"
#include <vector>
#include <map>
namespace Cantera
{
// Statistical mechanics
/*
* @ingroup spthermo
*/
//! Empty constructor
StatMech::StatMech()
: m_lowT(0.1), m_highT (1.0),
m_Pref(1.0E5), m_index (0) {}
// constructor used in templated instantiations
/*
* @param n Species index
* @param tlow Minimum temperature
* @param thigh Maximum temperature
* @param pref reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
StatMech::StatMech(int n, doublereal tlow, doublereal thigh,
doublereal pref,
const doublereal* coeffs,
std::string my_name) :
m_lowT (tlow),
m_highT (thigh),
m_Pref (pref),
m_index (n),
sp_name (my_name)
{
// should error on zero -- cannot take ln(0)
if(m_lowT <= 0.0){
throw CanteraError("Error in StatMech.cpp",
" Cannot take 0 tmin as input. \n\n");
}
buildmap();
}
// copy constructor
/*
* @param b object to be copied
*/
StatMech::StatMech(const StatMech& b) :
m_lowT (b.m_lowT),
m_highT (b.m_highT),
m_Pref (b.m_Pref),
m_index (b.m_index)
{
}
// assignment operator
/*
* @param b object to be copied
*/
StatMech& StatMech::operator=(const StatMech& b) {
if (&b != this) {
m_lowT = b.m_lowT;
m_highT = b.m_highT;
m_Pref = b.m_Pref;
m_index = b.m_index;
}
return *this;
}
// Destructor
StatMech::~StatMech() {
}
// duplicator
SpeciesThermoInterpType *
StatMech::duplMyselfAsSpeciesThermoInterpType() const {
StatMech* np = new StatMech(*this);
return (SpeciesThermoInterpType *) np;
}
// Returns the minimum temperature that the thermo
// parameterization is valid
doublereal StatMech::minTemp() const
{
return m_lowT;
}
// Returns the maximum temperature that the thermo
// parameterization is valid
doublereal StatMech::maxTemp() const {
return m_highT;
}
// Returns the reference pressure (Pa)
doublereal StatMech::refPressure() const { return m_Pref; }
// Returns an integer representing the type of parameterization
int StatMech::reportType() const {
return STAT;
}
// Returns an integer representing the species index
size_t StatMech::speciesIndex() const {
return m_index;
}
int StatMech::buildmap()
{
// build vector of strings
std::vector<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
int ii;
for(ii=0; ii < SS.size(); ii++)
{
name_map[SS[ii]]=(new species);
// init to crazy defaults
name_map[SS[ii]]->nvib = -1;
name_map[SS[ii]]->cfs = -1;
name_map[SS[ii]]->mol_weight = -1;
name_map[SS[ii]]->theta[0] =0.0;
name_map[SS[ii]]->theta[1] =0.0;
name_map[SS[ii]]->theta[2] =0.0;
name_map[SS[ii]]->theta[3] =0.0;
name_map[SS[ii]]->theta[4] =0.0;
}
// now set all species information
// build Air
name_map["Air"]->cfs = 2.5;
name_map["Air"]->mol_weight=28.96;
name_map["Air"]->nvib=0;
// build CPAir
name_map["CPAir"]->cfs = 2.5;
name_map["CPAir"]->mol_weight=28.96;
name_map["CPAir"]->nvib=0;
// build Ar
name_map["Ar"]->cfs = 1.5;
name_map["Ar"]->mol_weight=39.944;
name_map["Ar"]->nvib=0;
// build Ar+
name_map["Ar+"]->cfs = 1.5;
name_map["Ar+"]->mol_weight=39.94345;
name_map["Ar+"]->nvib=0;
// build C
name_map["C"]->cfs = 1.5;
name_map["C"]->mol_weight=12.011;
name_map["C"]->nvib=0;
// build C+
name_map["C+"]->cfs = 1.5;
name_map["C+"]->mol_weight=12.01045;
name_map["C+"]->nvib=0;
// C2
name_map["C2"]->cfs=2.5;
name_map["C2"]->mol_weight=24.022;
name_map["C2"]->nvib=1;
name_map["C2"]->theta[0]=2.6687e3;
// C2H
name_map["C2H"]->cfs=2.5;
name_map["C2H"]->mol_weight=25.03;
name_map["C2H"]->nvib=3;
name_map["C2H"]->theta[0]=5.20100e+03;
name_map["C2H"]->theta[1]=7.20000e+03;
name_map["C2H"]->theta[2]=2.66100e+03;
// C2H2
name_map["C2H2"]->cfs=2.5;
name_map["C2H2"]->mol_weight=26.038;
name_map["C2H2"]->nvib=5;
name_map["C2H2"]->theta[0]=4.85290e+03;
name_map["C2H2"]->theta[1]=2.84000e+03;
name_map["C2H2"]->theta[2]=4.72490e+03;
name_map["C2H2"]->theta[3]=8.81830e+02;
name_map["C2H2"]->theta[4]=1.05080e+03;
// C3
name_map["C3"]->cfs=2.5;
name_map["C3"]->mol_weight=36.033;
name_map["C3"]->nvib=3;
name_map["C3"]->theta[0]=1.84500e+03;
name_map["C3"]->theta[1]=7.78700e+02;
name_map["C3"]->theta[2]=3.11760e+03;
// CF
name_map["CF"]->cfs=2.5;
name_map["CF"]->mol_weight=31.00940;
name_map["CF"]->nvib=1;
name_map["CF"]->theta[0]=1.88214e+03;
// CF2
name_map["CF2"]->cfs=3;
name_map["CF2"]->mol_weight=50.00780;
name_map["CF2"]->nvib=3;
name_map["CF2"]->theta[0]=1.76120e+03;
name_map["CF2"]->theta[1]=9.56820e+02;
name_map["CF2"]->theta[2]=1.60000e+03;
// CF3
name_map["CF3"]->cfs=3;
name_map["CF3"]->mol_weight=69.00620;
name_map["CF3"]->nvib=4;
name_map["CF3"]->theta[0]=1.56800e+03;
name_map["CF3"]->theta[1]=1.00900e+03;
name_map["CF3"]->theta[2]=1.81150e+03;
name_map["CF3"]->theta[3]=7.36680e+02;
// CF4
name_map["CF4"]->cfs=3;
name_map["CF4"]->mol_weight=88.00460;
name_map["CF4"]->nvib=4;
name_map["CF4"]->theta[0]=1.30720e+03;
name_map["CF4"]->theta[1]=6.25892e+02;
name_map["CF4"]->theta[2]=1.84540e+03;
name_map["CF4"]->theta[3]=9.08950e+02;
// CH
name_map["CH"]->cfs=2.5;
name_map["CH"]->mol_weight=13.01900;
name_map["CH"]->nvib=1;
name_map["CH"]->theta[0]=4.11290e+03;
// CH2
name_map["CH2"]->cfs=3;
name_map["CH2"]->mol_weight=14.02700;
name_map["CH2"]->nvib=3;
name_map["CH2"]->theta[0]=4.31650e+03;
name_map["CH2"]->theta[1]=1.95972e+03;
name_map["CH2"]->theta[2]=4.60432e+03;
// CH3
name_map["CH3"]->cfs=3;
name_map["CH3"]->mol_weight=15.03500;
name_map["CH3"]->nvib=4;
name_map["CH3"]->theta[0]=4.31650e+03;
name_map["CH3"]->theta[1]=8.73370e+02;
name_map["CH3"]->theta[2]=4.54960e+03;
name_map["CH3"]->theta[3]=2.01150e+03;
// CH4
name_map["CH4"]->cfs=3;
name_map["CH4"]->mol_weight=16.04300;
name_map["CH4"]->nvib=4;
name_map["CH4"]->theta[0]=4.19660e+03;
name_map["CH4"]->theta[1]=2.20620e+03;
name_map["CH4"]->theta[2]=4.34450e+03;
name_map["CH4"]->theta[3]=1.88600e+03;
// Cl
name_map["Cl"]->cfs=1.5;
name_map["Cl"]->mol_weight=35.45300;
name_map["Cl"]->nvib=0;
// Cl2
name_map["Cl2"]->cfs=2.5;
name_map["Cl2"]->mol_weight=70.96;
name_map["Cl2"]->nvib=1;
name_map["Cl2"]->theta[0]=8.05355e+02;
// CN
name_map["CN"]->cfs=2.5;
name_map["CN"]->mol_weight=26.01900;
name_map["CN"]->nvib=1;
name_map["CN"]->theta[0]=2.97610e+03;
// CN+
name_map["CN+"]->cfs=2.5;
name_map["CN+"]->mol_weight=26.01845;
name_map["CN+"]->nvib=1;
name_map["CN+"]->theta[0]=2.92520e+03;
// CO
name_map["CO"]->cfs=2.5;
name_map["CO"]->mol_weight=28.01100;
name_map["CO"]->nvib=1;
name_map["CO"]->theta[0]=3.12200e+03;
// CO+
name_map["CO+"]->cfs=2.5;
name_map["CO+"]->mol_weight=28.01045;
name_map["CO+"]->nvib=1;
name_map["CO+"]->theta[0]=3.18800e+03;
// CO2
name_map["CO2"]->cfs=2.5;
name_map["CO2"]->mol_weight=44.01100;
name_map["CO2"]->nvib=3;
name_map["CO2"]->theta[0]=1.91870e+03;
name_map["CO2"]->theta[1]=9.59660e+02;
name_map["CO2"]->theta[2]=3.38210e+03;
// F
name_map["F"]->cfs=1.5;
name_map["F"]->mol_weight=18.99840;
name_map["F"]->nvib=0;
// F2
name_map["F2"]->cfs=2.5;
name_map["F2"]->mol_weight=37.99680;
name_map["F2"]->nvib=1;
name_map["F2"]->theta[0]=1.32020e+03;
// H
name_map["H"]->cfs=1.5;
name_map["H"]->mol_weight=1;
name_map["H"]->nvib=0;
// H+
name_map["H+"]->cfs=1.5;
name_map["H+"]->mol_weight=1.00745;
name_map["H+"]->nvib=0;
// H2
name_map["H2"]->cfs=2.5;
name_map["H2"]->mol_weight=2.01600;
name_map["H2"]->nvib=1;
name_map["H2"]->theta[0]=6.33140e+03;
// H2+
name_map["H2+"]->cfs=2.5;
name_map["H2+"]->mol_weight=2.01545;
name_map["H2+"]->nvib=1;
name_map["H2+"]->theta[0]=3.34280e+03;
// H2O
name_map["H2O"]->cfs=3.0;
name_map["H2O"]->mol_weight=18.01600;
name_map["H2O"]->nvib=3;
name_map["H2O"]->theta[0]=5.26130e+03;
name_map["H2O"]->theta[1]=2.29460e+03;
name_map["H2O"]->theta[2]=5.40395e+03;
// HCl
name_map["HCl"]->cfs=2.5;
name_map["HCl"]->mol_weight=36.46100;
name_map["HCl"]->nvib=1;
name_map["HCl"]->theta[0]=4.30330e+03;
// HCN
name_map["HCN"]->cfs=2.5;
name_map["HCN"]->mol_weight=27.02700;
name_map["HCN"]->nvib=3;
name_map["HCN"]->theta[0]=3.01620e+03;
name_map["HCN"]->theta[1]=1.02660e+03;
name_map["HCN"]->theta[2]=4.76450e+03;
// He
name_map["He"]->cfs=1.5;
name_map["He"]->mol_weight=4.00300;
name_map["He"]->nvib=0;
// He+
name_map["He+"]->cfs=1.5;
name_map["He+"]->mol_weight=4.00245;
name_map["He+"]->nvib=0;
// N
name_map["N"]->cfs=1.5;
name_map["N"]->mol_weight=14.008;
name_map["N"]->nvib=0;
// Ne
name_map["Ne"]->cfs=1.5;
name_map["Ne"]->mol_weight=20.17900;
name_map["Ne"]->nvib=0;
// N+
name_map["N+"]->cfs=1.5;
name_map["N+"]->mol_weight=14.00745;
name_map["N+"]->nvib=0;
// N2
name_map["N2"]->cfs=2.5;
name_map["N2"]->mol_weight=28.01600;
name_map["N2"]->nvib=1;
name_map["N2"]->theta[0]=3.39500e+03;
// N2+
name_map["N2+"]->cfs=2.5;
name_map["N2+"]->mol_weight=28.01545;
name_map["N2+"]->nvib=1;
name_map["N2+"]->theta[0]=3.17580e+03;
// CPN2
name_map["CPN2"]->cfs=2.5;
name_map["CPN2"]->mol_weight=28.01600;
name_map["CPN2"]->nvib=0;
// NCO
name_map["NCO"]->cfs=2.5;
name_map["NCO"]->mol_weight=42.01900;
name_map["NCO"]->nvib=3;
name_map["NCO"]->theta[0]=1.83600e+03;
name_map["NCO"]->theta[1]=7.67100e+02;
name_map["NCO"]->theta[2]=2.76800e+03;
// NH
name_map["NH"]->cfs=2.5;
name_map["NH"]->mol_weight=15.01600;
name_map["NH"]->nvib=1;
name_map["NH"]->theta[0]=4.72240e+03;
// NH+
name_map["NH+"]->cfs=2.5;
name_map["NH+"]->mol_weight=15.01545;
name_map["NH+"]->nvib=0;
// NH2
name_map["NH2"]->cfs=2.5;
name_map["NH2"]->mol_weight=16.02400;
name_map["NH2"]->nvib=0;
// NH3
name_map["NH3"]->cfs=2.5;
name_map["NH3"]->mol_weight=17.03200;
name_map["NH3"]->nvib=4;
name_map["NH3"]->theta[0]=4.78100e+03;
name_map["NH3"]->theta[1]=1.47040e+03;
name_map["NH3"]->theta[2]=4.95440e+03;
name_map["NH3"]->theta[3]=2.34070e+03;
// NO
name_map["NO"]->cfs=2.5;
name_map["NO"]->mol_weight=30.00800;
name_map["NO"]->nvib=1;
name_map["NO"]->theta[0]=2.81700e+03;
// NO+
name_map["NO+"]->cfs=2.5;
name_map["NO+"]->mol_weight=30.00745;
name_map["NO+"]->nvib=1;
name_map["NO+"]->theta[0]=3.42100e+03;
// NO2
name_map["NO2"]->cfs=3;
name_map["NO2"]->mol_weight=46.00800;
name_map["NO2"]->nvib=3;
name_map["NO2"]->theta[0]=1.07900e+03;
name_map["NO2"]->theta[1]=1.90000e+03;
name_map["NO2"]->theta[2]=2.32700e+03;
// O
name_map["O"]->cfs=1.5;
name_map["O"]->mol_weight=16.000;
name_map["O"]->nvib=0;
// O+
name_map["O+"]->cfs=1.5;
name_map["O+"]->mol_weight=15.99945;
name_map["O+"]->nvib=0;
// O2
name_map["O2"]->cfs=2.5;
name_map["O2"]->mol_weight=32.00000;
name_map["O2"]->nvib=1;
name_map["O2"]->theta[0]=2.23900e+03;
// O2
name_map["O2+"]->cfs=2.5;
name_map["O2+"]->mol_weight=31.99945;
name_map["O2+"]->nvib=1;
name_map["O2+"]->theta[0]=2.74120e+03;
// OH
name_map["OH"]->cfs=2.5;
name_map["OH"]->mol_weight=17.00800;
name_map["OH"]->nvib=1;
name_map["OH"]->theta[0]=5.37820e+03;
// Si
name_map["Si"]->cfs=1.5;
name_map["Si"]->mol_weight=28.08550;
name_map["Si"]->nvib=0;
// SiO
name_map["SiO"]->cfs=2.5;
name_map["SiO"]->mol_weight=44.08550;
name_map["SiO"]->nvib=1;
name_map["SiO"]->theta[0]=1.78640e+03;
// electron
name_map["e"]->cfs=1.5;
name_map["e"]->mol_weight=0.00055;
name_map["e"]->nvib=0;
int dum = 0;
for(ii=0; ii < SS.size(); ii++)
{
// check nvib was initalized for all species
if(name_map[SS[ii]]->nvib == -1)
{
std::cout << name_map[SS[ii]]->nvib << std::endl;
throw CanteraError("Error in StatMech.cpp",
"nvib not initialized!. \n\n");
}
else
{
// check that theta is initalized
for(int i=0;i<name_map[SS[ii]]->nvib;i++)
{
if(name_map[SS[ii]]->theta[i] <= 0.0)
{
throw CanteraError("Error in StatMech.cpp",
"theta not initalized!. \n\n");
}
}
// check that no non-zero theta exist
// for any theta larger than nvib!
for(int i=name_map[SS[ii]]->nvib;i<5;i++)
{
if(name_map[SS[ii]]->theta[i] != 0.0)
{
std::string err = "bad theta value for "+SS[ii]+"\n";
throw CanteraError("StatMech.cpp",err);
}
} // done with for loop
}
// check mol weight was initialized for all species
if(name_map[SS[ii]]->mol_weight == -1)
{
std::cout << name_map[SS[ii]]->mol_weight << std::endl;
throw CanteraError("Error in StatMech.cpp",
"mol_weight not initialized!. \n\n");
}
// cfs was initialized for all species
if(name_map[SS[ii]]->cfs == -1)
{
std::cout << name_map[SS[ii]]->cfs << std::endl;
throw CanteraError("Error in StatMech.cpp",
"cfs not initialized!. \n\n");
}
} // done with sanity checks
// mark it zero, dude
return 0;
}
// Update the properties for this species
/**
*
* \f[
* \frac{C_p^0(T)}{R} = \frac{C_v^0(T)}{R} + 1
* \f]
*
* Where,
* \f[
* \frac{C_v^0(T)}{R} = \frac{C_v^{tr}(T)}{R} + \frac{C_v^{vib}(T)}{R}
* \f]
*
*
* @param tt vector of temperature polynomials
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
void StatMech::updateProperties(const doublereal* tt,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const {
std::map<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;
}
// Compute the reference-state property of one species
/*
* Given temperature T in K, this method updates the values of
* the non-dimensional heat capacity at constant pressure,
* enthalpy, and entropy, at the reference pressure, Pref
* of one of the species. The species index is used
* to reference into the cp_R, h_RT, and s_R arrays.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param temp Temperature (Kelvin)
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
void StatMech::updatePropertiesTemp(const doublereal temp,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const {
double tPoly[1];
tPoly[0] = temp;
updateProperties(tPoly, cp_R, h_RT, s_R);
}
//This utility function reports back the type of
// parameterization and all of the parameters for the
// species, index.
/*
* All parameters are output variables
*
* @param n Species index
* @param type Integer type of the standard type
* @param tlow output - Minimum temperature
* @param thigh output - Maximum temperature
* @param pref output - reference pressure (Pa).
* @param coeffs Vector of species state data
*/
void StatMech::reportParameters(size_t &n, int &type,
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const
{
species* s;
n = m_index;
type = STAT;
tlow = m_lowT;
thigh = m_highT;
pref = m_Pref;
for (int i = 0; i < 9; i++)
{
coeffs[i] = 0.0;
}
doublereal temp = coeffs[0];
coeffs[1] = m_lowT;
coeffs[2] = m_highT;
// get species name, to gather species properties
// pointer to map location of particular species
if(name_map.find(sp_name) != name_map.end())
{
s = name_map.find(sp_name)->second;
}
else
{
//std::cout << sp_name << std::endl;
throw CanteraError("StatMech.cpp",
"species properties not found!. \n\n");
}
double theta = 0.0;
doublereal cvib = 0;
// vibrational energy
for(int i=0; i< s->nvib; i++)
{
theta = s->theta[i];
cvib += GasConstant * theta * (theta* exp(theta/temp)/(temp*temp))/((exp(theta/temp)-1) * (exp(theta/temp)-1));
}
// load vibrational energy
coeffs[3] = GasConstant * s->cfs;
coeffs[4] = cvib;
}
// Modify parameters for the standard state
/*
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
void StatMech::modifyParameters(doublereal* coeffs)
{
}
}

4
src/transport/Makefile.am
View file

@ -6,8 +6,8 @@ cc_sources = AqueousTransport.cpp LiquidTransport.cpp MMCollisionInt.cpp \
LiquidTranInteraction.cpp LiquidTransportData.cpp \
LiquidTransportParams.cpp TortuosityBase.cpp \
TortuosityBruggeman.cpp TortuosityMaxwell.cpp \
TortuosityPercolation.cpp TransportParams.cpp GasTransport.cpp
TortuosityPercolation.cpp TransportParams.cpp \
GasTransport.cpp PecosTransport.cpp
AM_CPPFLAGS = -I$(top_builddir)/include
AM_CXXFLAGS = $(AM_CPPFLAGS)

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

@ -0,0 +1,758 @@
/**
* @file PecosTransport.cpp
* Mixture-averaged transport properties.
*
*/
/* $Author$
* $Revision$
* $Date$
*/
#include "cantera/thermo/ThermoPhase.h"
#include "cantera/transport/PecosTransport.h"
#include "cantera/base/utilities.h"
#include "cantera/transport/TransportParams.h"
#include "cantera/transport/TransportFactory.h"
#include "cantera/base/stringUtils.h"
#include "cantera/thermo/IdealGasPhase.h"
#include <iostream>
using namespace std;
/**
* Mole fractions below MIN_X will be set to MIN_X when computing
* transport properties.
*/
#define MIN_X 1.e-20
namespace Cantera {
//////////////////// class PecosTransport methods //////////////
PecosTransport::PecosTransport() :
m_nsp(0),
m_tmin(-1.0),
m_tmax(100000.),
m_temp(-1.0),
m_logt(0.0)
{
}
bool PecosTransport::initGas( GasTransportParams& tr ) {
// constant substance attributes
m_thermo = tr.thermo;
m_nsp = m_thermo->nSpecies();
m_tmin = m_thermo->minTemp();
m_tmax = m_thermo->maxTemp();
// make a local copy of the molecular weights
m_mw.resize(m_nsp);
copy(m_thermo->molecularWeights().begin(),
m_thermo->molecularWeights().end(), m_mw.begin());
// copy polynomials and parameters into local storage
m_poly = tr.poly;
m_visccoeffs = tr.visccoeffs;
m_condcoeffs = tr.condcoeffs;
m_diffcoeffs = tr.diffcoeffs;
m_zrot = tr.zrot;
m_crot = tr.crot;
m_epsilon = tr.epsilon;
m_mode = tr.mode_;
m_diam = tr.diam;
m_eps = tr.eps;
m_alpha = tr.alpha;
m_dipoleDiag.resize(m_nsp);
for (int i = 0; i < m_nsp; i++) {
m_dipoleDiag[i] = tr.dipole(i,i);
}
m_phi.resize(m_nsp, m_nsp, 0.0);
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
int j, k;
for (j = 0; j < m_nsp; j++)
for (k = j; k < m_nsp; k++) {
m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
m_wratjk(k,j) = sqrt(m_wratjk(j,k));
m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
}
m_polytempvec.resize(5);
m_visc.resize(m_nsp);
m_sqvisc.resize(m_nsp);
m_cond.resize(m_nsp);
m_bdiff.resize(m_nsp, m_nsp);
m_molefracs.resize(m_nsp);
m_spwork.resize(m_nsp);
// set flags all false
m_viscmix_ok = false;
m_viscwt_ok = false;
m_spvisc_ok = false;
m_spcond_ok = false;
m_condmix_ok = false;
m_spcond_ok = false;
m_diffmix_ok = false;
m_abc_ok = false;
// read blottner fit parameters (A,B,C)
cout << "reading blottner";
read_blottner_transport_table();
cout << "done with blottner";
// set specific heats
cv_rot.resize(m_nsp);
cp_R.resize(m_nsp);
cv_int.resize(m_nsp);
for (k = 0; k < m_nsp; k++) {
cv_rot[k] = tr.crot[k];
cp_R[k] = ((IdealGasPhase*)tr.thermo)->cp_R_ref()[k];
cv_int[k] = cp_R[k] - 2.5 - cv_rot[k];
}
return true;
}
/*********************************************************
*
* Public methods
*
*********************************************************/
/****************** viscosity ******************************/
/**
* The viscosity is computed using the Wilke mixture rule.
* \f[
* \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}.
* \f]
* Here \f$ \mu_k \f$ is the viscosity of pure species \e k,
* and
* \f[
* \Phi_{k,j} = \frac{\left[1
* + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2}
* {\sqrt{8}\sqrt{1 + M_k/M_j}}
* \f]
* @see updateViscosity_T();
*/
doublereal PecosTransport::viscosity() {
update_T();
update_C();
if (m_viscmix_ok) return m_viscmix;
doublereal vismix = 0.0;
int k;
// update m_visc and m_phi if necessary
if (!m_viscwt_ok) updateViscosity_T();
multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork));
for (k = 0; k < m_nsp; k++) {
vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
}
m_viscmix = vismix;
return vismix;
}
/******************* binary diffusion coefficients **************/
/*
*
* Using Ramshaw's self-consistent Effective Binary Diffusion
* (1990, J. Non-Equilib. Thermo)
* Adding more doxygen would be good here
*/
void PecosTransport::getBinaryDiffCoeffs(const int ld, doublereal* const d) {
int i,j;
update_T();
// if necessary, evaluate the binary diffusion coefficents
if (!m_bindiff_ok) updateDiff_T();
doublereal rp = 1.0/pressure_ig();
for (i = 0; i < m_nsp; i++)
for (j = 0; j < m_nsp; j++) {
d[ld*j + i] = rp * m_bdiff(i,j);
}
}
void PecosTransport::getMobilities(doublereal* const mobil) {
int k;
getMixDiffCoeffs(DATA_PTR(m_spwork));
doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
for (k = 0; k < m_nsp; k++) {
mobil[k] = c1 * m_spwork[k] * m_thermo->charge(k);
}
}
/****************** thermal conductivity **********************/
/**
* The thermal conductivity is computed using the Wilke mixture rule.
* \f[
* \k = \sum_s \frac{k_s X_s}{\sum_j \Phi_{s,j} X_j}.
* \f]
* Here \f$ \k_s \f$ is the conductivity of pure species \e s,
* and
* \f[
* \Phi_{s,j} = \frac{\left[1
* + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_s}}\right)}\right]^2}
* {\sqrt{8}\sqrt{1 + M_s/M_j}}
* \f]
* @see updateCond_T();
*/
doublereal PecosTransport::thermalConductivity() {
int k;
doublereal lambda = 0.0;
update_T();
update_C();
// update m_cond and m_phi if necessary
if (!m_spcond_ok) updateCond_T();
if (!m_condmix_ok) {
multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork));
for (k = 0; k < m_nsp; k++) {
lambda += m_molefracs[k] * m_cond[k]/m_spwork[k]; //denom;
}
}
m_lambda = lambda;
return m_lambda;
}
/****************** thermal diffusion coefficients ************/
/**
* Thermal diffusion is not considered in this pecos
* model. To include thermal diffusion, use transport manager
* MultiTransport instead. This methods fills out array dt with
* zeros.
*/
void PecosTransport::getThermalDiffCoeffs(doublereal* const dt) {
int k;
for (k = 0; k < m_nsp; k++) {
dt[k] = 0.0;
}
}
/**
* @param ndim The number of spatial dimensions (1, 2, or 3).
* @param grad_T The temperature gradient (ignored in this model).
* @param ldx Leading dimension of the grad_X array.
* The diffusive mass flux of species \e k is computed from
* \f[
* \vec{j}_k = -n M_k D_k \nabla X_k + \frac{\rho_k}{\rho} \sum_r n M_r D_r \nabla X_r
* \f]
*
* This is neglective pressure, forced and thermal diffusion.
*
*/
void PecosTransport::getSpeciesFluxes(int ndim,
const doublereal* grad_T, int ldx, const doublereal* grad_X,
int ldf, doublereal* fluxes) {
int n, k;
update_T();
update_C();
getMixDiffCoeffs(DATA_PTR(m_spwork));
const vector_fp& mw = m_thermo->molecularWeights();
const doublereal* y = m_thermo->massFractions();
doublereal rhon = m_thermo->molarDensity();
vector_fp sum(ndim,0.0);
doublereal correction=0.0;
// grab 2nd (summation) term -- still need to multiply by mass fraction (\rho_s / \rho)
for (k = 0; k < m_nsp; k++)
{
correction += rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k];
}
for (n = 0; n < ndim; n++) {
for (k = 0; k < m_nsp; k++) {
fluxes[n*ldf + k] = -rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k] + y[k]*correction;
sum[n] += fluxes[n*ldf + k];
}
}
// add correction flux to enforce sum to zero
for (n = 0; n < ndim; n++) {
for (k = 0; k < m_nsp; k++) {
fluxes[n*ldf + k] -= y[k]*sum[n];
}
}
}
/**
* Mixture-averaged diffusion coefficients [m^2/s].
*
* For the single species case or the pure fluid case
* the routine returns the self-diffusion coefficient.
* This is need to avoid a Nan result in the formula
* below.
*/
void PecosTransport::getMixDiffCoeffs(doublereal* const d) {
update_T();
update_C();
// update the binary diffusion coefficients if necessary
if (!m_bindiff_ok) updateDiff_T();
int k, j;
doublereal mmw = m_thermo->meanMolecularWeight();
doublereal sumxw = 0.0, sum2;
doublereal p = pressure_ig();
if (m_nsp == 1) {
d[0] = m_bdiff(0,0) / p;
} else {
for (k = 0; k < m_nsp; k++) sumxw += m_molefracs[k] * m_mw[k];
for (k = 0; k < m_nsp; k++) {
sum2 = 0.0;
for (j = 0; j < m_nsp; j++) {
if (j != k) {
sum2 += m_molefracs[j] / m_bdiff(j,k);
}
}
if (sum2 <= 0.0) {
d[k] = m_bdiff(k,k) / p;
} else {
d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
}
}
}
}
void PecosTransport::getMixDiffCoeffsMole(doublereal* const d)
{
update_T();
update_C();
// update the binary diffusion coefficients if necessary
if (!m_bindiff_ok) {
updateDiff_T();
}
doublereal p = m_thermo->pressure();
if (m_nsp == 1) {
d[0] = m_bdiff(0,0) / p;
} else {
for (size_t k = 0; k < m_nsp; k++) {
double sum2 = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != k) {
sum2 += m_molefracs[j] / m_bdiff(j,k);
}
}
if (sum2 <= 0.0) {
d[k] = m_bdiff(k,k) / p;
} else {
d[k] = (1 - m_molefracs[k]) / (p * sum2);
}
}
}
}
void PecosTransport::getMixDiffCoeffsMass(doublereal* const d)
{
update_T();
update_C();
// update the binary diffusion coefficients if necessary
if (!m_bindiff_ok) {
updateDiff_T();
}
doublereal mmw = m_thermo->meanMolecularWeight();
doublereal p = m_thermo->pressure();
if (m_nsp == 1) {
d[0] = m_bdiff(0,0) / p;
} else {
for (size_t k=0; k<m_nsp; k++) {
double sum1 = 0.0;
double sum2 = 0.0;
for (size_t 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;
}
/**
* @internal This is called the first time any transport property
* is requested from Mixture after the concentrations
* have changed.
*/
void PecosTransport::update_C()
{
// signal that concentration-dependent quantities will need to
// be recomputed before use, and update the local mole
// fractions.
m_viscmix_ok = false;
m_diffmix_ok = false;
m_condmix_ok = false;
m_thermo->getMoleFractions(DATA_PTR(m_molefracs));
// add an offset to avoid a pure species condition
int k;
for (k = 0; k < m_nsp; k++) {
m_molefracs[k] = std::max(MIN_X, m_molefracs[k]);
}
}
/*************************************************************************
*
* methods to update temperature-dependent properties
*
*************************************************************************/
/**
*
* Update the temperature-dependent parts of the mixture-averaged
* thermal conductivity.
*
* Calculated as,
* \f[
* k= \mu_s (5/2 * C_{v,s}^{trans} + C_{v,s}^{rot} + C_{v,s}^{vib}
* \f]
*
*
*/
void PecosTransport::updateCond_T() {
int k;
doublereal fivehalves = 5/2;
for (k = 0; k < m_nsp; k++) {
// need to add cv_elec in the future
m_cond[k] = m_visc[k] * ( fivehalves * cv_int[k] + cv_rot[k] + m_thermo->cv_vib(k,m_temp) );
}
m_spcond_ok = true;
m_condmix_ok = false;
}
/**
* Update the binary diffusion coefficients. These are evaluated
* from the polynomial fits at unit pressure (1 Pa).
*/
void PecosTransport::updateDiff_T() {
// evaluate binary diffusion coefficients at unit pressure
int i,j;
int ic = 0;
if (m_mode == CK_Mode) {
for (i = 0; i < m_nsp; i++) {
for (j = i; j < m_nsp; j++) {
m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic]));
m_bdiff(j,i) = m_bdiff(i,j);
ic++;
}
}
}
else {
for (i = 0; i < m_nsp; i++) {
for (j = i; j < m_nsp; j++) {
m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec,
m_diffcoeffs[ic]);
m_bdiff(j,i) = m_bdiff(i,j);
ic++;
}
}
}
m_bindiff_ok = true;
m_diffmix_ok = false;
}
/**
*
* Update the pure-species viscosities. (Pa-s) = (kg/m/sec)
*
* Using Blottner fit for viscosity. Defines kinematic viscosity
* of the form
* \f[
* \mu_s\left(T\right) = 0.10 \exp\left(A_s\left(\log T\right)^2 + B_s\log T + C_s\right)
* \f]
* where \f$ A_s \f$, \f$ B_s \f$, and \f$ C_s \f$ are constants.
*
*/
void PecosTransport::updateSpeciesViscosities() {
// blottner
// return 0.10*std::exp(_a*(logT*logT) + _b*logT + _c);
int k;
// iterate over species, update pure-species viscosity
for (k = 0; k < m_nsp; k++) {
m_visc[k] = 0.10*std::exp(a[k]*(m_logt*m_logt) + b[k]*m_logt + c[k]);
m_sqvisc[k] = sqrt(m_visc[k]);
}
// time to update mixing
m_spvisc_ok = true;
}
/*
* read_blottner_transport_table()
* loads up A B and C for blottner fits
* hardcoded for air, will need to generalize later
*/
void PecosTransport::read_blottner_transport_table()
{
// istringstream blot
// ("Air 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
// "CPAir 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
// "N 1.15572000000e-02 6.03167900000e-01 -1.24327495000e+01\n"
// "N2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
// "CPN2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
// "NO 4.36378000000e-02 -3.35511000000e-02 -9.57674300000e+00\n"
// "O 2.03144000000e-02 4.29440400000e-01 -1.16031403000e+01\n"
// "O2 4.49290000000e-02 -8.26158000000e-02 -9.20194750000e+00\n"
// "C -8.3285e-3 0.7703240 -12.7378000\n"
// "C2 -8.4311e-3 0.7876060 -13.0268000\n"
// "C3 -8.4312e-3 0.7876090 -12.8240000\n"
// "C2H -2.4241e-2 1.0946550 -14.5835500\n"
// "CN -8.3811e-3 0.7860330 -12.9406000\n"
// "CO -0.019527394 1.013295 -13.97873\n"
// "CO2 -0.019527387 1.047818 -14.32212\n"
// "HCN -2.4241e-2 1.0946550 -14.5835500\n"
// "H -8.3912e-3 0.7743270 -13.6653000\n"
// "H2 -8.3346e-3 0.7815380 -13.5351000\n"
// "e 0.00000000000e+00 0.00000000000e+00 -1.16031403000e+01\n");
//
// from: AIAA-1997-2474 and Sandia Report SC-RR-70-754
//
// # Air -- Identical to N2 fit
// # N -- Sandia Report SC-RR-70-754
// # N2 -- Sandia Report SC-RR-70-754
// # CPN2 -- Identical to N2 fit
// # NO -- Sandia Report SC-RR-70-754
// # O -- Sandia Report SC-RR-70-754
// # O2 -- Sandia Report SC-RR-70-754
// # C -- AIAA-1997-2474
// # C2 -- AIAA-1997-2474
// # C3 -- AIAA-1997-2474
// # C2H -- wild-ass guess: identical to HCN fit
// # CN -- AIAA-1997-2474
// # CO -- AIAA-1997-2474
// # CO2 -- AIAA-1997-2474
// # HCN -- AIAA-1997-2474
// # H -- AIAA-1997-2474
// # H2 -- AIAA-1997-2474
// # e -- Sandia Report SC-RR-70-754
istringstream blot
("Air 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
"CPAir 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
"N 1.15572000000e-02 6.03167900000e-01 -1.24327495000e+01\n"
"N2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
"CPN2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
"NO 4.36378000000e-02 -3.35511000000e-02 -9.57674300000e+00\n"
"O 2.03144000000e-02 4.29440400000e-01 -1.16031403000e+01\n"
"O2 4.49290000000e-02 -8.26158000000e-02 -9.20194750000e+00\n"
"C -8.3285e-3 0.7703240 -12.7378000\n"
"C2 -8.4311e-3 0.7876060 -13.0268000\n"
"C3 -8.4312e-3 0.7876090 -12.8240000\n"
"C2H -2.4241e-2 1.0946550 -14.5835500\n"
"CN -8.3811e-3 0.7860330 -12.9406000\n"
"CO -0.019527394 1.013295 -13.97873\n"
"CO2 -0.019527387 1.047818 -14.32212\n"
"HCN -2.4241e-2 1.0946550 -14.5835500\n"
"H -8.3912e-3 0.7743270 -13.6653000\n"
"H2 -8.3346e-3 0.7815380 -13.5351000\n"
"e 0.00000000000e+00 0.00000000000e+00 -1.16031403000e+01\n");
string line;
string name;
string ss1,ss2,ss3,ss4,sss;
int k;
int i = 0;
while (std::getline(blot, line))
{
istringstream ss(line);
std::getline(ss, ss1, ' ');
std::getline(ss, ss2, ' ');
std::getline(ss, ss3, ' ');
std::getline(ss, ss4, ' ');
name = ss1;
// now put coefficients in correct species
for (k = 0; k < m_nsp; k++)
{
string sss = m_thermo->speciesName(k);
// this is the right species index
if(sss.compare(ss1) == 0)
{
a[k] = atof(ss2.c_str());
b[k] = atof(ss3.c_str());
c[k] = atof(ss4.c_str());
// index
i++;
}
else // default to air
{
a[k] = 0.026;
b[k] = 0.3;
c[k] = -11.3;
}
} // done with for loop
}
// for (k = 0; k < m_nsp; k++)
// {
// string sss = m_thermo->speciesName(k);
// cout << sss << endl;
// cout << a[k] << endl;
// cout << b[k] << endl;
// cout << c[k] << endl;
// }
// simple sanity check
// if(i != m_nsp-1)
// {
// std::cout << "error\n" << i << std::endl;
// }
}
/**
*
* Update the temperature-dependent viscosity terms.
* Updates the array of pure species viscosities, and the
* weighting functions in the viscosity mixture rule.
* The flag m_visc_ok is set to true.
*
*/
void PecosTransport::updateViscosity_T() {
doublereal vratiokj, wratiojk, factor1;
if (!m_spvisc_ok) updateSpeciesViscosities();
// see Eq. (9-5.15) of Reid, Prausnitz, and Poling
int j, k;
for (j = 0; j < m_nsp; j++) {
for (k = j; k < m_nsp; k++) {
vratiokj = m_visc[k]/m_visc[j];
wratiojk = m_mw[j]/m_mw[k];
// Note that m_wratjk(k,j) holds the square root of
// m_wratjk(j,k)!
factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
m_phi(k,j) = factor1*factor1 /
(SqrtEight * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
}
m_viscwt_ok = true;
}
// /**
// *
// * This function returns a Transport data object for a given species.
// *
// */
// struct GasTransportData PecosTransport::
// getGasTransportData(int kSpecies)
// {
// struct GasTransportData td;
// td.speciesName = m_thermo->speciesName(kSpecies);
// td.geometry = 2;
// if (m_crot[kSpecies] == 0.0) {
// td.geometry = 0;
// } else if (m_crot[kSpecies] == 1.0) {
// td.geometry = 1;
// }
// td.wellDepth = m_eps[kSpecies] / Boltzmann;
// td.dipoleMoment = m_dipoleDiag[kSpecies] * 1.0E25 / SqrtTen;
// td.diameter = m_diam(kSpecies, kSpecies) * 1.0E10;
// td.polarizability = m_alpha[kSpecies] * 1.0E30;
// td.rotRelaxNumber = m_zrot[kSpecies];
// return td;
// }
}

7
src/transport/TransportFactory.cpp
View file

@ -8,6 +8,7 @@
// known transport models
#include "cantera/transport/MultiTransport.h"
#include "cantera/transport/PecosTransport.h"
#include "cantera/transport/MixTransport.h"
#include "cantera/transport/SolidTransport.h"
#include "cantera/transport/DustyGasTransport.h"
@ -202,6 +203,7 @@ TransportFactory::TransportFactory() :
m_models["Aqueous"] = cAqueousTransport;
m_models["Simple"] = cSimpleTransport;
m_models["User"] = cUserTransport;
m_models["Pecos"] = cPecosTransport;
m_models["None"] = None;
//m_models["Radiative"] = cRadiative;
@ -371,6 +373,11 @@ Transport* TransportFactory::newTransport(std::string transportModel,
tr = new MixTransport;
initTransport(tr, phase, CK_Mode, log_level);
break;
// adding pecos transport model 2/13/12
case cPecosTransport:
tr = new PecosTransport;
initTransport(tr, phase, 0, log_level);
break;
case cSolidTransport:
tr = new SolidTransport;
tr->setThermo(*phase);

10
test_problems/PecosTransport/PecosTransport.cpp
View file

@ -32,14 +32,18 @@ using namespace std;
/*****************************************************************/
/*****************************************************************/
#include "cantera/Cantera.h"
#include "cantera/transport.h"
#include "cantera/IdealGasMix.h"
#include "cantera/transport/TransportFactory.h"
//#include "Cantera.h"
//#include "transport.h"
//#include "IdealGasMix.h"
//#include "TransportFactory.h"
using namespace Cantera;
using namespace Cantera_CXX;
//using namespace Cantera_CXX;
void printDbl(double val) {
if (fabs(val) < 5.0E-17) {