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

[Transport] Move gas transport initialization to class GasTransport

Browse source 
Model-specific initialization code belongs in the classes implementing
the corresponding model, not in TransportFactory.
This commit is contained in:
Ray Speth 2015-01-07 23:30:57 +00:00
parent 141e7d9b40
commit 07fd8e26c0
4 changed files with 786 additions and 812 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

150
include/cantera/transport/GasTransport.h
View file

@ -11,6 +11,8 @@
namespace Cantera
{
class MMCollisionInt;
//! Class GasTransport implements some functions and properties that are
//! shared by the MixTransport and MultiTransport classes.
//! @ingroup tranprops
@ -105,6 +107,19 @@ public:
*/
virtual void getMixDiffCoeffsMass(doublereal* const d);
//! Initialize a transport manager
/*!
* This routine sets up a gas-phase transport manager. It calculates the
* collision integrals and calls the initGas() function to populate the
* species-dependent data structure.
*
* @param thermo Pointer to the ThermoPhase object
* @param mode Chemkin compatible mode or not. This alters the
* specification of the collision integrals. defaults to no.
* @param log_level Defaults to zero, no logging
*/
virtual void init(thermo_t* thermo, int mode=0, int log_level=0);
protected:
GasTransport(ThermoPhase* thermo=0);
@ -133,10 +148,139 @@ protected:
//! Update the binary diffusion coefficients
/*!
* These are evaluated from the polynomial fits of the temperature at the unit pressure of 1 Pa.
* These are evaluated from the polynomial fits of the temperature at the
* unit pressure of 1 Pa.
*/
virtual void updateDiff_T();
//! @name Initialization
//! @{
//! Prepare to build a new kinetic-theory-based transport manager for
//! low-density gases
/*!
* This class fills up the GastransportParams structure for the current phase
*
* Uses polynomial fits to Monchick & Mason collision integrals. Store them
* in tr.
*
* @param transport_database Reference to a vector of pointers containing
* the transport database for each species
* @param thermo Pointer to the ThermoPhase object
* @param mode Mode -> Either it's CK_Mode, chemkin compatibility
* mode, or it is not We usually run with chemkin
* compatibility mode turned off.
* @param log_level log level
* @param tr GasTransportParams structure to be filled up with
* information
*/
void setupMM(const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int mode, int log_level,
GasTransportParams& tr);
//! Read the transport database
/*!
* Read transport property data from a file for a list of species. Given the
* name of a file containing transport property parameters and a list of
* species names, this method returns an instance of TransportParams
* containing the transport data for these species read from the file.
*
* @param thermo The phase with species corresponding to the transport data
* @param xspecies Vector of pointers to species XML_Node databases.
* @param log reference to an XML_Node that will contain the log (unused)
* @param names vector of species names that must be filled in with
* valid transport parameters
* @param tr Output object containing the transport parameters for
* the species listed in names (in the order of their
* listing in names).
*/
void getTransportData(const ThermoPhase& thermo,
const std::vector<const XML_Node*> &xspecies,
XML_Node& log, const std::vector<std::string>& names,
GasTransportParams& tr);
//! Corrections for polar-nonpolar binary diffusion coefficients
/*!
* Calculate corrections to the well depth parameter and the diameter for
* use in computing the binary diffusion coefficient of polar-nonpolar
* pairs. For more information about this correction, see Dixon-Lewis, Proc.
* Royal Society (1968).
*
* @param i Species one - this is a bimolecular correction routine
* @param j species two - this is a bimolecular correction routine
* @param tr Database of species properties read in from the input xml file.
* @param f_eps Multiplicative correction factor to be applied to epsilon(i,j)
* @param f_sigma Multiplicative correction factor to be applied to diam(i,j)
*/
void makePolarCorrections(size_t i, size_t j,
const GasTransportParams& tr, doublereal& f_eps,
doublereal& f_sigma);
//! Generate polynomial fits to collision integrals
/*!
* @param tr Reference to the GasTransportParams object that will
* contain the results.
* @param integrals interpolator for the collision integrals
*/
void fitCollisionIntegrals(GasTransportParams& tr,
MMCollisionInt& integrals);
//! Generate polynomial fits to the viscosity, conductivity, and
//! the binary diffusion coefficients
/*!
* If CK_mode, then the fits are of the form
* \f[
* \log(\eta(i)) = \sum_{n = 0}^3 a_n(i) (\log T)^n
* \f]
* and
* \f[
* \log(D(i,j)) = \sum_{n = 0}^3 a_n(i,j) (\log T)^n
* \f]
* Otherwise the fits are of the form
* \f[
* \eta(i)/sqrt(k_BT) = \sum_{n = 0}^4 a_n(i) (\log T)^n
* \f]
* and
* \f[
* D(i,j)/sqrt(k_BT)) = \sum_{n = 0}^4 a_n(i,j) (\log T)^n
* \f]
*
* @param tr Reference to the GasTransportParams object that will
* contain the results.
* @param integrals interpolator for the collision integrals
*/
void fitProperties(GasTransportParams& tr, MMCollisionInt& integrals);
//! Second-order correction to the binary diffusion coefficients
/*!
* Calculate second-order corrections to binary diffusion coefficient pair
* (dkj, djk). At first order, the binary diffusion coefficients are
* independent of composition, and d(k,j) = d(j,k). But at second order,
* there is a weak dependence on composition, with the result that d(k,j) !=
* d(j,k). This method computes the multiplier by which the first-order
* binary diffusion coefficient should be multiplied to produce the value
* correct to second order. The expressions here are taken from Marerro and
* Mason, J. Phys. Chem. Ref. Data, vol. 1, p. 3 (1972).
*
* @param t Temperature (K)
* @param tr Transport parameters
* @param integrals interpolator for the collision integrals
* @param k index of first species
* @param j index of second species
* @param xk Mole fraction of species k
* @param xj Mole fraction of species j
* @param fkj multiplier for d(k,j)
* @param fjk multiplier for d(j,k)
*
* @note This method is not used currently.
*/
void getBinDiffCorrection(doublereal t, const GasTransportParams& tr,
MMCollisionInt& integrals, size_t k,
size_t j, doublereal xk, doublereal xj,
doublereal& fkj, doublereal& fjk);
//! @}
//! Vector of species mole fractions. These are processed so that all mole
//! fractions are >= *Tiny*. Length = m_kk.
vector_fp m_molefracs;
@ -242,6 +386,10 @@ protected:
//! Matrix of binary diffusion coefficients at the reference pressure and
//! the current temperature Size is nsp x nsp.
DenseMatrix m_bdiff;
//! Boolean indicating whether to turn on verbose printing during
//! initialization
bool m_verbose;
};
} // namespace Cantera

136
include/cantera/transport/TransportFactory.h
View file

@ -17,9 +17,6 @@
namespace Cantera
{
// forward references
class MMCollisionInt;
//! Factory class for creating new instances of classes derived from Transport.
/*!
* Creates 'transport managers', which are classes derived from class
@ -110,20 +107,6 @@ public:
virtual Transport*
newTransport(thermo_t* thermo, int log_level=0);
//! Initialize an existing transport manager
/*!
* This routine sets up an existing gas-phase transport manager. It
* calculates the collision integrals and calls the initGas() function to
* populate the species-dependent data structure.
*
* @param tr Pointer to the Transport manager
* @param thermo Pointer to the ThermoPhase object
* @param mode Chemkin compatible mode or not. This alters the specification of the
* collision integrals. defaults to no.
* @param log_level Defaults to zero, no logging
*/
virtual void initTransport(Transport* tr, thermo_t* thermo, int mode=0, int log_level=0);
//! Initialize an existing transport manager for liquid phase
/*!
* This routine sets up an existing liquid-phase transport manager. It is
@ -166,26 +149,6 @@ private:
*/
TransportFactory();
//! Read the transport database
/*!
* Read transport property data from a file for a list of species.
* Given the name of a file containing transport property
* parameters and a list of species names, this method returns an
* instance of TransportParams containing the transport data for
* these species read from the file.
*
* @param thermo The phase with species corresponding to the transport data
* @param xspecies Vector of pointers to species XML_Node databases.
* @param log reference to an XML_Node that will contain the log (unused)
* @param names vector of species names that must be filled in with valid transport parameters
* @param tr Output object containing the transport parameters
* for the species listed in names (in the order of their listing
* in names).
*/
void getTransportData(const ThermoPhase& thermo, const std::vector<const XML_Node*> &xspecies,
XML_Node& log, const std::vector<std::string>& names,
GasTransportParams& tr);
//! Read transport property data from a file for a list of species that comprise
//! the phase.
/*!
@ -240,56 +203,6 @@ private:
const std::string phaseName,
SolidTransportData& tr);
//! Generate polynomial fits to the viscosity, conductivity, and
//! the binary diffusion coefficients
/*!
* If CK_mode, then the fits are of the form
* \f[
* \log(\eta(i)) = \sum_{n = 0}^3 a_n(i) (\log T)^n
* \f]
* and
* \f[
* \log(D(i,j)) = \sum_{n = 0}^3 a_n(i,j) (\log T)^n
* \f]
* Otherwise the fits are of the form
* \f[
* \eta(i)/sqrt(k_BT) = \sum_{n = 0}^4 a_n(i) (\log T)^n
* \f]
* and
* \f[
* D(i,j)/sqrt(k_BT)) = \sum_{n = 0}^4 a_n(i,j) (\log T)^n
* \f]
*
* @param tr Reference to the GasTransportParams object that will contain the results.
* @param integrals interpolator for the collision integrals
*/
void fitProperties(GasTransportParams& tr, MMCollisionInt& integrals);
//! Generate polynomial fits to collision integrals
/*!
* @param tr Reference to the GasTransportParams object that will contain the results.
* @param integrals interpolator for the collision integrals
*/
void fitCollisionIntegrals(GasTransportParams& tr,
MMCollisionInt& integrals);
//! Prepare to build a new kinetic-theory-based transport manager for low-density gases
/*!
* This class fills up the GastransportParams structure for the current phase
*
* Uses polynomial fits to Monchick & Mason collision integrals. store then in tr
*
* @param transport_database Reference to a vector of pointers containing the
* transport database for each species
* @param thermo Pointer to the ThermoPhase object
* @param mode Mode -> Either it's CK_Mode, chemkin compatibility mode, or it is not
* We usually run with chemkin compatibility mode turned off.
* @param log_level log level
* @param tr GasTransportParams structure to be filled up with information
*/
void setupMM(const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int mode, int log_level, GasTransportParams& tr);
//! Prepare to build a new transport manager for liquids assuming that
//! viscosity transport data is provided in Arrhenius form.
/*!
@ -307,55 +220,6 @@ private:
*/
void setupSolidTransport(thermo_t* thermo, int log_level, SolidTransportData& trParam);
//! Second-order correction to the binary diffusion coefficients
/*!
* Calculate second-order corrections to binary diffusion
* coefficient pair (dkj, djk). At first order, the binary
* diffusion coefficients are independent of composition, and
* d(k,j) = d(j,k). But at second order, there is a weak
* dependence on composition, with the result that d(k,j) !=
* d(j,k). This method computes the multiplier by which the
* first-order binary diffusion coefficient should be multiplied
* to produce the value correct to second order. The expressions
* here are taken from Marerro and Mason, J. Phys. Chem. Ref. Data, vol. 1, p. 3 (1972).
*
* @param t Temperature (K)
* @param tr Transport parameters
* @param integrals interpolator for the collision integrals
* @param k index of first species
* @param j index of second species
* @param xk Mole fraction of species k
* @param xj Mole fraction of species j
* @param fkj multiplier for d(k,j)
* @param fjk multiplier for d(j,k)
*
* @note This method is not used currently.
*/
void getBinDiffCorrection(doublereal t, const GasTransportParams& tr,
MMCollisionInt& integrals, size_t k,
size_t j, doublereal xk, doublereal xj,
doublereal& fkj, doublereal& fjk);
//! Corrections for polar-nonpolar binary diffusion coefficients
/*!
* Calculate corrections to the well depth parameter and the
* diameter for use in computing the binary diffusion coefficient
* of polar-nonpolar pairs. For more information about this
* correction, see Dixon-Lewis, Proc. Royal Society (1968).
*
* @param i Species one - this is a bimolecular correction routine
* @param j species two - this is a bimolecular correction routine
* @param tr Database of species properties read in from the input xml file.
* @param f_eps Multiplicative correction factor to be applied to epsilon(i,j)
* @param f_sigma Multiplicative correction factor to be applied to diam(i,j)
*/
void makePolarCorrections(size_t i, size_t j,
const GasTransportParams& tr, doublereal& f_eps,
doublereal& f_sigma);
//! Boolean indicating whether to turn on verbose printing
bool m_verbose;
//! Mapping between between the string name for a transport model and the
//! integer name.
std::map<std::string, int> m_models;

631
src/transport/GasTransport.cpp
View file

@ -1,10 +1,18 @@
//! @file GasTransport.cpp
#include "cantera/transport/GasTransport.h"
#include "cantera/transport/TransportParams.h"
#include "MMCollisionInt.h"
#include "cantera/base/ctml.h"
#include "cantera/base/stringUtils.h"
#include "cantera/numerics/polyfit.h"
namespace Cantera
{
//! polynomial degree used for fitting collision integrals
//! except in CK mode, where the degree is 6.
#define COLL_INT_POLY_DEGREE 8
GasTransport::GasTransport(ThermoPhase* thermo) :
Transport(thermo),
m_molefracs(0),
@ -31,7 +39,8 @@ GasTransport::GasTransport(ThermoPhase* thermo) :
m_t14(0.0),
m_t32(0.0),
m_diffcoeffs(0),
m_bdiff(0, 0)
m_bdiff(0, 0),
m_verbose(false)
{
}
@ -60,7 +69,8 @@ GasTransport::GasTransport(const GasTransport& right) :
m_t14(0.0),
m_t32(0.0),
m_diffcoeffs(0),
m_bdiff(0, 0)
m_bdiff(0, 0),
m_verbose(false)
{
}
@ -90,6 +100,7 @@ GasTransport& GasTransport::operator=(const GasTransport& right)
m_t32 = right.m_t32;
m_diffcoeffs = right.m_diffcoeffs;
m_bdiff = right.m_bdiff;
m_verbose = right.m_verbose;
return *this;
}
@ -372,4 +383,620 @@ void GasTransport::getMixDiffCoeffsMass(doublereal* const d)
}
}
void GasTransport::init(thermo_t* thermo, int mode, int log_level)
{
GasTransportParams trParam;
if (log_level == 0) {
m_verbose = 0;
}
// set up Monchick and Mason collision integrals
setupMM(thermo->speciesData(), thermo, mode, log_level, trParam);
// do model-specific initialization
initGas(trParam);
}
void GasTransport::setupMM(const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int mode, int log_level, GasTransportParams& tr)
{
// constant mixture attributes
tr.thermo = thermo;
tr.nsp_ = tr.thermo->nSpecies();
size_t nsp = tr.nsp_;
tr.tmin = thermo->minTemp();
tr.tmax = thermo->maxTemp();
tr.mw.resize(nsp);
tr.log_level = log_level;
copy(tr.thermo->molecularWeights().begin(),
tr.thermo->molecularWeights().end(), tr.mw.begin());
tr.mode_ = mode;
tr.epsilon.resize(nsp, nsp, 0.0);
tr.delta.resize(nsp, nsp, 0.0);
tr.reducedMass.resize(nsp, nsp, 0.0);
tr.dipole.resize(nsp, nsp, 0.0);
tr.diam.resize(nsp, nsp, 0.0);
tr.crot.resize(nsp);
tr.zrot.resize(nsp);
tr.polar.resize(nsp, false);
tr.alpha.resize(nsp, 0.0);
tr.poly.resize(nsp);
tr.sigma.resize(nsp);
tr.eps.resize(nsp);
tr.w_ac.resize(nsp);
XML_Node root, log;
getTransportData(*thermo, transport_database, log, tr.thermo->speciesNames(), tr);
for (size_t i = 0; i < nsp; i++) {
tr.poly[i].resize(nsp);
}
double tstar_min = 1.e8, tstar_max = 0.0;
double f_eps, f_sigma;
for (size_t i = 0; i < nsp; i++) {
for (size_t j = i; j < nsp; j++) {
// the reduced mass
tr.reducedMass(i,j) = tr.mw[i] * tr.mw[j] / (Avogadro * (tr.mw[i] + tr.mw[j]));
// hard-sphere diameter for (i,j) collisions
tr.diam(i,j) = 0.5*(tr.sigma[i] + tr.sigma[j]);
// the effective well depth for (i,j) collisions
tr.epsilon(i,j) = sqrt(tr.eps[i]*tr.eps[j]);
// The polynomial fits of collision integrals vs. T*
// will be done for the T* from tstar_min to tstar_max
tstar_min = std::min(tstar_min, Boltzmann * tr.tmin/tr.epsilon(i,j));
tstar_max = std::max(tstar_max, Boltzmann * tr.tmax/tr.epsilon(i,j));
// the effective dipole moment for (i,j) collisions
tr.dipole(i,j) = sqrt(tr.dipole(i,i)*tr.dipole(j,j));
// reduced dipole moment delta* (nondimensional)
double d = tr.diam(i,j);
tr.delta(i,j) = 0.5 * tr.dipole(i,j)*tr.dipole(i,j)
/ (4 * Pi * epsilon_0 * tr.epsilon(i,j) * d * d * d);
makePolarCorrections(i, j, tr, f_eps, f_sigma);
tr.diam(i,j) *= f_sigma;
tr.epsilon(i,j) *= f_eps;
// properties are symmetric
tr.reducedMass(j,i) = tr.reducedMass(i,j);
tr.diam(j,i) = tr.diam(i,j);
tr.epsilon(j,i) = tr.epsilon(i,j);
tr.dipole(j,i) = tr.dipole(i,j);
tr.delta(j,i) = tr.delta(i,j);
}
}
// Chemkin fits the entire T* range in the Monchick and Mason tables,
// so modify tstar_min and tstar_max if in Chemkin compatibility mode
if (mode == CK_Mode) {
tstar_min = 0.101;
tstar_max = 99.9;
}
// initialize the collision integral calculator for the desired T* range
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** collision_integrals ***\n");
}
MMCollisionInt integrals;
integrals.init(tstar_min, tstar_max, log_level);
fitCollisionIntegrals(tr, integrals);
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** end of collision_integrals ***\n");
}
// make polynomial fits
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** property fits ***\n");
}
fitProperties(tr, integrals);
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** end of property fits ***\n");
}
}
void GasTransport::getTransportData(const ThermoPhase& thermo,
const std::vector<const XML_Node*> &xspecies,
XML_Node& log,
const std::vector<std::string> &names,
GasTransportParams& tr)
{
std::map<std::string, size_t> speciesIndices;
for (size_t i = 0; i < names.size(); i++) {
speciesIndices[names[i]] = i;
}
for (size_t i = 0; i < xspecies.size(); i++) {
const XML_Node& sp = *xspecies[i];
// Find the index for this species in 'names'
size_t j = getValue(speciesIndices, sp["name"], npos);
if (j == npos) {
// Don't need transport data for this species
continue;
}
XML_Node& node = sp.child("transport");
// parameters are converted to SI units before storing
double nAtoms = 0;
size_t kSpec = thermo.speciesIndex(sp["name"]);
for (size_t m = 0; m < thermo.nElements(); m++) {
nAtoms += thermo.nAtoms(kSpec, m);
}
// Molecular geometry; rotational heat capacity / R
XML_Node* geomNode = ctml::getByTitle(node, "geometry");
std::string geom = (geomNode) ? geomNode->value() : "";
if (geom == "atom") {
if (nAtoms != 1) {
throw CanteraError("GasTransport::getTransportData",
"invalid geometry. 'atom' specified,"
" but species contains multiple atoms.");
}
tr.crot[j] = 0.0;
} else if (geom == "linear") {
if (nAtoms == 1) {
throw CanteraError("GasTransport::getTransportData",
"invalid geometry. 'linear' specified,"
" but species only contains one atom.");
}
tr.crot[j] = 1.0;
} else if (geom == "nonlinear") {
if (nAtoms < 3) {
throw CanteraError("GasTransport::getTransportData",
"invalid geometry. 'nonlinear' specified,"
" but species only contains " + fp2str(nAtoms) + " atoms.");
}
tr.crot[j] = 1.5;
} else {
throw CanteraError("GasTransport::getTransportData",
"invalid geometry");
}
// Pitzer's acentric factor:
double acentric;
ctml::getOptionalFloat(node, "acentric_factor", acentric);
if (acentric) {
tr.w_ac[j] = acentric;
}
// Well-depth parameter in Kelvin (converted to Joules)
double welldepth = ctml::getFloat(node, "LJ_welldepth");
if (welldepth >= 0.0) {
tr.eps[j] = Boltzmann * welldepth;
} else {
throw CanteraError("GasTransport::getTransportData",
"negative well depth");
}
// Lennard-Jones diameter of the molecule, given in Angstroms.
double diam = ctml::getFloat(node, "LJ_diameter");
if (diam > 0.0) {
tr.sigma[j] = 1.e-10 * diam; // A -> m
} else {
throw CanteraError("GasTransport::getTransportData",
"negative or zero diameter");
}
// Dipole moment of the molecule.
// Given in Debye (a Debye is 1e-18 statC-m or 3.3356e-30 C-m)
double dipole = ctml::getFloat(node, "dipoleMoment");
if (dipole >= 0.0) {
tr.dipole(j,j) = 1e-21 / lightSpeed * dipole;
tr.polar[j] = (dipole > 0.0);
} else {
throw CanteraError("GasTransport::getTransportData",
"negative dipole moment");
}
// Polarizability of the molecule, given in cubic Angstroms.
double polar = ctml::getFloat(node, "polarizability");
if (polar >= 0.0) {
tr.alpha[j] = 1.e-30 * polar; // A^3 -> m^3
} else {
throw CanteraError("GasTransport::getTransportData",
"negative polarizability");
}
// Rotational relaxation number. (Number of collisions it takes to
// equilibrate the rotational dofs with the temperature)
double rot = ctml::getFloat(node, "rotRelax");
if (rot >= 0.0) {
tr.zrot[j] = std::max(1.0, rot);
} else {
throw CanteraError("GasTransport::getTransportData",
"negative rotation relaxation number");
}
}
}
void GasTransport::makePolarCorrections(size_t i, size_t j,
const GasTransportParams& tr, doublereal& f_eps, doublereal& f_sigma)
{
// no correction if both are nonpolar, or both are polar
if (tr.polar[i] == tr.polar[j]) {
f_eps = 1.0;
f_sigma = 1.0;
return;
}
// corrections to the effective diameter and well depth
// if one is polar and one is non-polar
size_t kp = (tr.polar[i] ? i : j); // the polar one
size_t knp = (i == kp ? j : i); // the nonpolar one
double d3np, d3p, alpha_star, mu_p_star, xi;
d3np = pow(tr.sigma[knp],3);
d3p = pow(tr.sigma[kp],3);
alpha_star = tr.alpha[knp]/d3np;
mu_p_star = tr.dipole(kp,kp)/sqrt(4 * Pi * epsilon_0 * d3p * tr.eps[kp]);
xi = 1.0 + 0.25 * alpha_star * mu_p_star * mu_p_star *
sqrt(tr.eps[kp]/tr.eps[knp]);
f_sigma = pow(xi, -1.0/6.0);
f_eps = xi*xi;
}
void GasTransport::fitCollisionIntegrals(GasTransportParams& tr,
MMCollisionInt& integrals)
{
vector_fp::iterator dptr;
double dstar;
size_t nsp = tr.nsp_;
int mode = tr.mode_;
// Chemkin fits to sixth order polynomials
int degree = (mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("tstar_fits\n"
"fits to A*, B*, and C* vs. log(T*).\n"
"These are done only for the required dstar(j,k) values.\n\n");
if (tr.log_level < 3) {
writelog("*** polynomial coefficients not printed (log_level < 3) ***\n");
}
}
for (size_t i = 0; i < nsp; i++) {
for (size_t j = i; j < nsp; j++) {
// Chemkin fits only delta* = 0
if (mode != CK_Mode) {
dstar = tr.delta(i,j);
} else {
dstar = 0.0;
}
// if a fit has already been generated for delta* = tr.delta(i,j),
// then use it. Otherwise, make a new fit, and add tr.delta(i,j) to
// the list of delta* values for which fits have been done.
// 'find' returns a pointer to end() if not found
dptr = find(tr.fitlist.begin(), tr.fitlist.end(), dstar);
if (dptr == tr.fitlist.end()) {
vector_fp ca(degree+1), cb(degree+1), cc(degree+1);
vector_fp co22(degree+1);
integrals.fit(degree, dstar,
DATA_PTR(ca), DATA_PTR(cb), DATA_PTR(cc));
integrals.fit_omega22(degree, dstar,
DATA_PTR(co22));
tr.omega22_poly.push_back(co22);
tr.astar_poly.push_back(ca);
tr.bstar_poly.push_back(cb);
tr.cstar_poly.push_back(cc);
tr.poly[i][j] = static_cast<int>(tr.astar_poly.size()) - 1;
tr.fitlist.push_back(dstar);
}
// delta* found in fitlist, so just point to this polynomial
else {
tr.poly[i][j] = static_cast<int>((dptr - tr.fitlist.begin()));
}
tr.poly[j][i] = tr.poly[i][j];
}
}
}
void GasTransport::fitProperties(GasTransportParams& tr,
MMCollisionInt& integrals)
{
int ndeg = 0;
// number of points to use in generating fit data
const size_t np = 50;
int mode = tr.mode_;
int degree = (mode == CK_Mode ? 3 : 4);
double dt = (tr.tmax - tr.tmin)/(np-1);
vector_fp tlog(np), spvisc(np), spcond(np);
vector_fp w(np), w2(np);
// generate array of log(t) values
for (size_t n = 0; n < np; n++) {
double t = tr.tmin + dt*n;
tlog[n] = log(t);
}
// vector of polynomial coefficients
vector_fp c(degree + 1), c2(degree + 1);
// fit the pure-species viscosity and thermal conductivity for each species
if (DEBUG_MODE_ENABLED && tr.log_level < 2 && m_verbose) {
writelog("*** polynomial coefficients not printed (log_level < 2) ***\n");
}
double sqrt_T, visc, err, relerr,
mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0;
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("Polynomial fits for viscosity:\n");
if (mode == CK_Mode) {
writelog("log(viscosity) fit to cubic polynomial in log(T)\n");
} else {
writelogf("viscosity/sqrt(T) fit to polynomial of degree "
"%d in log(T)", degree);
}
}
double cp_R, cond, w_RT, f_int, A_factor, B_factor, c1, cv_rot, cv_int,
f_rot, f_trans, om11, diffcoeff;
for (size_t k = 0; k < tr.nsp_; k++) {
for (size_t n = 0; n < np; n++) {
double t = tr.tmin + dt*n;
tr.thermo->setTemperature(t);
vector_fp cp_R_all(tr.thermo->nSpecies());
tr.thermo->getCp_R_ref(&cp_R_all[0]);
cp_R = cp_R_all[k];
double tstar = Boltzmann * t/ tr.eps[k];
sqrt_T = sqrt(t);
double om22 = integrals.omega22(tstar, tr.delta(k,k));
om11 = integrals.omega11(tstar, tr.delta(k,k));
// self-diffusion coefficient, without polar corrections
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/tr.reducedMass(k,k)) *
pow((Boltzmann * t), 1.5)/
(Pi * tr.sigma[k] * tr.sigma[k] * om11);
// viscosity
visc = FiveSixteenths
* sqrt(Pi * tr.mw[k] * Boltzmann * t / Avogadro) /
(om22 * Pi * tr.sigma[k]*tr.sigma[k]);
// thermal conductivity
w_RT = tr.mw[k]/(GasConstant * t);
f_int = w_RT * diffcoeff/visc;
cv_rot = tr.crot[k];
A_factor = 2.5 - f_int;
B_factor = tr.zrot[k] + 2.0/Pi * (5.0/3.0 * cv_rot + f_int);
c1 = 2.0/Pi * A_factor/B_factor;
cv_int = cp_R - 2.5 - cv_rot;
f_rot = f_int * (1.0 + c1);
f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5);
cond = (visc/tr.mw[k])*GasConstant*(f_trans * 1.5
+ f_rot * cv_rot + f_int * cv_int);
if (mode == CK_Mode) {
spvisc[n] = log(visc);
spcond[n] = log(cond);
w[n] = -1.0;
w2[n] = -1.0;
} else {
// the viscosity should be proportional approximately to
// sqrt(T); therefore, visc/sqrt(T) should have only a weak
// temperature dependence. And since the mixture rule requires
// the square root of the pure-species viscosity, fit the square
// root of (visc/sqrt(T)) to avoid having to compute square
// roots in the mixture rule.
spvisc[n] = sqrt(visc/sqrt_T);
// the pure-species conductivity scales approximately with
// sqrt(T). Unlike the viscosity, there is no reason here to fit
// the square root, since a different mixture rule is used.
spcond[n] = cond/sqrt_T;
w[n] = 1.0/(spvisc[n]*spvisc[n]);
w2[n] = 1.0/(spcond[n]*spcond[n]);
}
}
polyfit(np, DATA_PTR(tlog), DATA_PTR(spvisc),
DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c));
polyfit(np, DATA_PTR(tlog), DATA_PTR(spcond),
DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c2));
// evaluate max fit errors for viscosity
for (size_t n = 0; n < np; n++) {
double val, fit;
if (mode == CK_Mode) {
val = exp(spvisc[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c)));
} else {
sqrt_T = exp(0.5*tlog[n]);
val = sqrt_T * pow(spvisc[n],2);
fit = sqrt_T * pow(poly4(tlog[n], DATA_PTR(c)),2);
}
err = fit - val;
relerr = err/val;
mxerr = std::max(mxerr, fabs(err));
mxrelerr = std::max(mxrelerr, fabs(relerr));
}
// evaluate max fit errors for conductivity
for (size_t n = 0; n < np; n++) {
double val, fit;
if (mode == CK_Mode) {
val = exp(spcond[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c2)));
} else {
sqrt_T = exp(0.5*tlog[n]);
val = sqrt_T * spcond[n];
fit = sqrt_T * poly4(tlog[n], DATA_PTR(c2));
}
err = fit - val;
relerr = err/val;
mxerr_cond = std::max(mxerr_cond, fabs(err));
mxrelerr_cond = std::max(mxrelerr_cond, fabs(relerr));
}
tr.visccoeffs.push_back(c);
tr.condcoeffs.push_back(c2);
if (DEBUG_MODE_ENABLED && tr.log_level >= 2 && m_verbose) {
writelog(tr.thermo->speciesName(k) + ": [" + vec2str(c) + "]\n");
}
}
if (DEBUG_MODE_ENABLED && m_verbose) {
writelogf("Maximum viscosity absolute error: %12.6g\n", mxerr);
writelogf("Maximum viscosity relative error: %12.6g\n", mxrelerr);
writelog("\nPolynomial fits for conductivity:\n");
if (mode == CK_Mode)
writelog("log(conductivity) fit to cubic polynomial in log(T)");
else {
writelogf("conductivity/sqrt(T) fit to "
"polynomial of degree %d in log(T)", degree);
}
if (tr.log_level >= 2)
for (size_t k = 0; k < tr.nsp_; k++) {
writelog(tr.thermo->speciesName(k) + ": [" +
vec2str(tr.condcoeffs[k]) + "]\n");
}
writelogf("Maximum conductivity absolute error: %12.6g\n", mxerr_cond);
writelogf("Maximum conductivity relative error: %12.6g\n", mxrelerr_cond);
// fit the binary diffusion coefficients for each species pair
writelogf("\nbinary diffusion coefficients:\n");
if (mode == CK_Mode)
writelog("log(D) fit to cubic polynomial in log(T)");
else {
writelogf("D/T**(3/2) fit to polynomial of degree %d in log(T)",degree);
}
}
mxerr = 0.0, mxrelerr = 0.0;
vector_fp diff(np + 1);
double eps, sigma;
for (size_t k = 0; k < tr.nsp_; k++) {
for (size_t j = k; j < tr.nsp_; j++) {
for (size_t n = 0; n < np; n++) {
double t = tr.tmin + dt*n;
eps = tr.epsilon(j,k);
double tstar = Boltzmann * t/eps;
sigma = tr.diam(j,k);
om11 = integrals.omega11(tstar, tr.delta(j,k));
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/tr.reducedMass(k,j)) *
pow(Boltzmann * t, 1.5) /
(Pi * sigma * sigma * om11);
// 2nd order correction
// NOTE: THIS CORRECTION IS NOT APPLIED
double fkj, fjk;
getBinDiffCorrection(t, tr, integrals, k, j, 1.0, 1.0, fkj, fjk);
if (mode == CK_Mode) {
diff[n] = log(diffcoeff);
w[n] = -1.0;
} else {
diff[n] = diffcoeff/pow(t, 1.5);
w[n] = 1.0/(diff[n]*diff[n]);
}
}
polyfit(np, DATA_PTR(tlog), DATA_PTR(diff),
DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c));
for (size_t n = 0; n < np; n++) {
double val, fit;
if (mode == CK_Mode) {
val = exp(diff[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c)));
} else {
double t = exp(tlog[n]);
double pre = pow(t, 1.5);
val = pre * diff[n];
fit = pre * poly4(tlog[n], DATA_PTR(c));
}
err = fit - val;
relerr = err/val;
mxerr = std::max(mxerr, fabs(err));
mxrelerr = std::max(mxrelerr, fabs(relerr));
}
tr.diffcoeffs.push_back(c);
if (DEBUG_MODE_ENABLED && tr.log_level >= 2 && m_verbose) {
writelog(tr.thermo->speciesName(k) + "__" +
tr.thermo->speciesName(j) + ": [" + vec2str(c) + "]\n");
}
}
}
if (DEBUG_MODE_ENABLED && m_verbose) {
writelogf("Maximum binary diffusion coefficient absolute error:"
" %12.6g\n", mxerr);
writelogf("Maximum binary diffusion coefficient relative error:"
"%12.6g", mxrelerr);
}
}
void GasTransport::getBinDiffCorrection(double t,
const GasTransportParams& tr, MMCollisionInt& integrals,
size_t k, size_t j, double xk, double xj, double& fkj, double& fjk)
{
double w1 = tr.mw[k];
double w2 = tr.mw[j];
double wsum = w1 + w2;
double wmwp = (w1 - w2)/wsum;
double sqw12 = sqrt(w1*w2);
double sig1 = tr.sigma[k];
double sig2 = tr.sigma[j];
double sig12 = 0.5*(tr.sigma[k] + tr.sigma[j]);
double sigratio = sig1*sig1/(sig2*sig2);
double sigratio2 = sig1*sig1/(sig12*sig12);
double sigratio3 = sig2*sig2/(sig12*sig12);
double tstar1 = Boltzmann * t / tr.eps[k];
double tstar2 = Boltzmann * t / tr.eps[j];
double tstar12 = Boltzmann * t / sqrt(tr.eps[k] * tr.eps[j]);
double om22_1 = integrals.omega22(tstar1, tr.delta(k,k));
double om22_2 = integrals.omega22(tstar2, tr.delta(j,j));
double om11_12 = integrals.omega11(tstar12, tr.delta(k,j));
double astar_12 = integrals.astar(tstar12, tr.delta(k,j));
double bstar_12 = integrals.bstar(tstar12, tr.delta(k,j));
double cstar_12 = integrals.cstar(tstar12, tr.delta(k,j));
double cnst = sigratio * sqrt(2.0*w2/wsum) * 2.0 * w1*w1/(wsum * w2);
double p1 = cnst * om22_1 / om11_12;
cnst = (1.0/sigratio) * sqrt(2.0*w1/wsum) * 2.0*w2*w2/(wsum*w1);
double p2 = cnst * om22_2 / om11_12;
double p12 = 15.0 * wmwp*wmwp + 8.0*w1*w2*astar_12/(wsum*wsum);
cnst = (2.0/(w2*wsum))*sqrt(2.0*w2/wsum)*sigratio2;
double q1 = cnst*((2.5 - 1.2*bstar_12)*w1*w1 + 3.0*w2*w2
+ 1.6*w1*w2*astar_12);
cnst = (2.0/(w1*wsum))*sqrt(2.0*w1/wsum)*sigratio3;
double q2 = cnst*((2.5 - 1.2*bstar_12)*w2*w2 + 3.0*w1*w1
+ 1.6*w1*w2*astar_12);
double q12 = wmwp*wmwp*15.0*(2.5 - 1.2*bstar_12)
+ 4.0*w1*w2*astar_12*(11.0 - 2.4*bstar_12)/(wsum*wsum)
+ 1.6*wsum*om22_1*om22_2/(om11_12*om11_12*sqw12)
* sigratio2 * sigratio3;
cnst = 6.0*cstar_12 - 5.0;
fkj = 1.0 + 0.1*cnst*cnst *
(p1*xk*xk + p2*xj*xj + p12*xk*xj)/
(q1*xk*xk + q2*xj*xj + q12*xk*xj);
fjk = 1.0 + 0.1*cnst*cnst *
(p2*xk*xk + p1*xj*xj + p12*xk*xj)/
(q2*xk*xk + q1*xj*xj + q12*xk*xj);
}
}

681
src/transport/TransportFactory.cpp
View file

@ -17,35 +17,19 @@
#include "cantera/transport/TransportFactory.h"
#include "cantera/transport/SolidTransportData.h"
#include "cantera/numerics/polyfit.h"
#include "MMCollisionInt.h"
#include "cantera/base/ctml.h"
#include "cantera/base/stringUtils.h"
#include "cantera/base/utilities.h"
using namespace std;
//! polynomial degree used for fitting collision integrals
//! except in CK mode, where the degree is 6.
#define COLL_INT_POLY_DEGREE 8
namespace Cantera
{
/////////////////////////// constants //////////////////////////
//@ \cond
const doublereal ThreeSixteenths = 3.0/16.0;
const doublereal TwoOverPi = 2.0/Pi;
const doublereal FiveThirds = 5.0/3.0;
//@ \endcond
TransportFactory* TransportFactory::s_factory = 0;
// declaration of static storage for the mutex
mutex_t TransportFactory::transport_mutex;
////////////////////////// exceptions /////////////////////////
//! Exception thrown if an error is encountered while reading the transport database
class TransportDBError : public CanteraError
{
@ -62,100 +46,7 @@ public:
//////////////////// class TransportFactory methods //////////////
void TransportFactory::getBinDiffCorrection(doublereal t,
const GasTransportParams& tr, MMCollisionInt& integrals,
size_t k, size_t j, doublereal xk, doublereal xj,
doublereal& fkj, doublereal& fjk)
{
doublereal w1, w2, wsum, sig1, sig2, sig12, sigratio, sigratio2,
sigratio3, tstar1, tstar2, tstar12,
om22_1, om22_2, om11_12, astar_12, bstar_12, cstar_12,
cnst, wmwp, sqw12, p1, p2, p12, q1, q2, q12;
w1 = tr.mw[k];
w2 = tr.mw[j];
wsum = w1 + w2;
wmwp = (w1 - w2)/wsum;
sqw12 = sqrt(w1*w2);
sig1 = tr.sigma[k];
sig2 = tr.sigma[j];
sig12 = 0.5*(tr.sigma[k] + tr.sigma[j]);
sigratio = sig1*sig1/(sig2*sig2);
sigratio2 = sig1*sig1/(sig12*sig12);
sigratio3 = sig2*sig2/(sig12*sig12);
tstar1 = Boltzmann * t / tr.eps[k];
tstar2 = Boltzmann * t / tr.eps[j];
tstar12 = Boltzmann * t / sqrt(tr.eps[k] * tr.eps[j]);
om22_1 = integrals.omega22(tstar1, tr.delta(k,k));
om22_2 = integrals.omega22(tstar2, tr.delta(j,j));
om11_12 = integrals.omega11(tstar12, tr.delta(k,j));
astar_12 = integrals.astar(tstar12, tr.delta(k,j));
bstar_12 = integrals.bstar(tstar12, tr.delta(k,j));
cstar_12 = integrals.cstar(tstar12, tr.delta(k,j));
cnst = sigratio * sqrt(2.0*w2/wsum) * 2.0 *
w1*w1/(wsum * w2);
p1 = cnst * om22_1 / om11_12;
cnst = (1.0/sigratio) * sqrt(2.0*w1/wsum) * 2.0*w2*w2/(wsum*w1);
p2 = cnst * om22_2 / om11_12;
p12 = 15.0 * wmwp*wmwp + 8.0*w1*w2*astar_12/(wsum*wsum);
cnst = (2.0/(w2*wsum))*sqrt(2.0*w2/wsum)*sigratio2;
q1 = cnst*((2.5 - 1.2*bstar_12)*w1*w1 + 3.0*w2*w2
+ 1.6*w1*w2*astar_12);
cnst = (2.0/(w1*wsum))*sqrt(2.0*w1/wsum)*sigratio3;
q2 = cnst*((2.5 - 1.2*bstar_12)*w2*w2 + 3.0*w1*w1
+ 1.6*w1*w2*astar_12);
q12 = wmwp*wmwp*15.0*(2.5 - 1.2*bstar_12)
+ 4.0*w1*w2*astar_12*(11.0 - 2.4*bstar_12)/(wsum*wsum)
+ 1.6*wsum*om22_1*om22_2/(om11_12*om11_12*sqw12)
* sigratio2 * sigratio3;
cnst = 6.0*cstar_12 - 5.0;
fkj = 1.0 + 0.1*cnst*cnst *
(p1*xk*xk + p2*xj*xj + p12*xk*xj)/
(q1*xk*xk + q2*xj*xj + q12*xk*xj);
fjk = 1.0 + 0.1*cnst*cnst *
(p2*xk*xk + p1*xj*xj + p12*xk*xj)/
(q2*xk*xk + q1*xj*xj + q12*xk*xj);
}
void TransportFactory::makePolarCorrections(size_t i, size_t j,
const GasTransportParams& tr, doublereal& f_eps, doublereal& f_sigma)
{
// no correction if both are nonpolar, or both are polar
if (tr.polar[i] == tr.polar[j]) {
f_eps = 1.0;
f_sigma = 1.0;
return;
}
// corrections to the effective diameter and well depth
// if one is polar and one is non-polar
size_t kp = (tr.polar[i] ? i : j); // the polar one
size_t knp = (i == kp ? j : i); // the nonpolar one
doublereal d3np, d3p, alpha_star, mu_p_star, xi;
d3np = pow(tr.sigma[knp],3);
d3p = pow(tr.sigma[kp],3);
alpha_star = tr.alpha[knp]/d3np;
mu_p_star = tr.dipole(kp,kp)/sqrt(4 * Pi * epsilon_0 * d3p * tr.eps[kp]);
xi = 1.0 + 0.25 * alpha_star * mu_p_star * mu_p_star *
sqrt(tr.eps[kp]/tr.eps[knp]);
f_sigma = pow(xi, -1.0/6.0);
f_eps = xi*xi;
}
TransportFactory::TransportFactory() :
m_verbose(false)
TransportFactory::TransportFactory()
{
m_models["Mix"] = cMixtureAveraged;
m_models["Multi"] = cMulticomponent;
@ -323,39 +214,34 @@ Transport* TransportFactory::newTransport(const std::string& transportModel,
break;
case cMulticomponent:
tr = new MultiTransport;
initTransport(tr, phase, 0, log_level);
dynamic_cast<GasTransport*>(tr)->init(phase, 0, log_level);
break;
case CK_Multicomponent:
tr = new MultiTransport;
initTransport(tr, phase, CK_Mode, log_level);
dynamic_cast<GasTransport*>(tr)->init(phase, CK_Mode, log_level);
break;
case cMixtureAveraged:
tr = new MixTransport;
initTransport(tr, phase, 0, log_level);
dynamic_cast<GasTransport*>(tr)->init(phase, 0, log_level);
break;
case CK_MixtureAveraged:
tr = new MixTransport;
initTransport(tr, phase, CK_Mode, log_level);
dynamic_cast<GasTransport*>(tr)->init(phase, CK_Mode, log_level);
break;
case cHighP:
tr = new HighPressureGasTransport;
initTransport(tr, phase, 0, log_level);
break;
// adding pecos transport model 2/13/12
case cPecosTransport:
tr = new PecosTransport;
initTransport(tr, phase, 0, log_level);
dynamic_cast<GasTransport*>(tr)->init(phase, 0, log_level);
break;
case cSolidTransport:
tr = new SolidTransport;
initSolidTransport(tr, phase, log_level);
tr->setThermo(*phase);
dynamic_cast<GasTransport*>(tr)->setThermo(*phase);
break;
case cDustyGasTransport:
tr = new DustyGasTransport;
gastr = new MultiTransport;
initTransport(gastr, phase, 0, log_level);
dynamic_cast<GasTransport*>(gastr)->init(phase, 0, log_level);
dtr = (DustyGasTransport*)tr;
dtr->initialize(phase, gastr);
break;
@ -400,117 +286,6 @@ Transport* TransportFactory::newTransport(thermo_t* phase, int log_level)
return newTransport(transportModel, phase,log_level);
}
void TransportFactory::setupMM(const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int mode, int log_level, GasTransportParams& tr)
{
// constant mixture attributes
tr.thermo = thermo;
tr.nsp_ = tr.thermo->nSpecies();
size_t nsp = tr.nsp_;
tr.tmin = thermo->minTemp();
tr.tmax = thermo->maxTemp();
tr.mw.resize(nsp);
tr.log_level = log_level;
copy(tr.thermo->molecularWeights().begin(), tr.thermo->molecularWeights().end(), tr.mw.begin());
tr.mode_ = mode;
tr.epsilon.resize(nsp, nsp, 0.0);
tr.delta.resize(nsp, nsp, 0.0);
tr.reducedMass.resize(nsp, nsp, 0.0);
tr.dipole.resize(nsp, nsp, 0.0);
tr.diam.resize(nsp, nsp, 0.0);
tr.crot.resize(nsp);
tr.zrot.resize(nsp);
tr.polar.resize(nsp, false);
tr.alpha.resize(nsp, 0.0);
tr.poly.resize(nsp);
tr.sigma.resize(nsp);
tr.eps.resize(nsp);
tr.w_ac.resize(nsp);
XML_Node root, log;
getTransportData(*thermo, transport_database, log, tr.thermo->speciesNames(), tr);
for (size_t i = 0; i < nsp; i++) {
tr.poly[i].resize(nsp);
}
doublereal tstar_min = 1.e8, tstar_max = 0.0;
doublereal f_eps, f_sigma;
DenseMatrix& diam = tr.diam;
DenseMatrix& epsilon = tr.epsilon;
for (size_t i = 0; i < nsp; i++) {
for (size_t j = i; j < nsp; j++) {
// the reduced mass
tr.reducedMass(i,j) = tr.mw[i] * tr.mw[j] / (Avogadro * (tr.mw[i] + tr.mw[j]));
// hard-sphere diameter for (i,j) collisions
diam(i,j) = 0.5*(tr.sigma[i] + tr.sigma[j]);
// the effective well depth for (i,j) collisions
epsilon(i,j) = sqrt(tr.eps[i]*tr.eps[j]);
// The polynomial fits of collision integrals vs. T*
// will be done for the T* from tstar_min to tstar_max
tstar_min = std::min(tstar_min, Boltzmann * tr.tmin/epsilon(i,j));
tstar_max = std::max(tstar_max, Boltzmann * tr.tmax/epsilon(i,j));
// the effective dipole moment for (i,j) collisions
tr.dipole(i,j) = sqrt(tr.dipole(i,i)*tr.dipole(j,j));
// reduced dipole moment delta* (nondimensional)
doublereal d = diam(i,j);
tr.delta(i,j) = 0.5 * tr.dipole(i,j)*tr.dipole(i,j)
/ (4 * Pi * epsilon_0 * epsilon(i,j) * d * d * d);
makePolarCorrections(i, j, tr, f_eps, f_sigma);
tr.diam(i,j) *= f_sigma;
epsilon(i,j) *= f_eps;
// properties are symmetric
tr.reducedMass(j,i) = tr.reducedMass(i,j);
diam(j,i) = diam(i,j);
epsilon(j,i) = epsilon(i,j);
tr.dipole(j,i) = tr.dipole(i,j);
tr.delta(j,i) = tr.delta(i,j);
}
}
// Chemkin fits the entire T* range in the Monchick and Mason tables,
// so modify tstar_min and tstar_max if in Chemkin compatibility mode
if (mode == CK_Mode) {
tstar_min = 0.101;
tstar_max = 99.9;
}
// initialize the collision integral calculator for the desired
// T* range
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** collision_integrals ***\n");
}
MMCollisionInt integrals;
integrals.init(tstar_min, tstar_max, log_level);
fitCollisionIntegrals(tr, integrals);
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** end of collision_integrals ***\n");
}
// make polynomial fits
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** property fits ***\n");
}
fitProperties(tr, integrals);
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("*** end of property fits ***\n");
}
}
void TransportFactory::setupLiquidTransport(thermo_t* thermo, int log_level,
LiquidTransportParams& trParam)
{
@ -585,23 +360,6 @@ void TransportFactory::setupSolidTransport(thermo_t* thermo, int log_level,
}
}
void TransportFactory::initTransport(Transport* tran,
thermo_t* thermo, int mode, int log_level)
{
ScopedLock transportLock(transport_mutex);
const std::vector<const XML_Node*> & transport_database = thermo->speciesData();
GasTransportParams trParam;
if (log_level == 0) {
m_verbose = 0;
}
// set up Monchick and Mason collision integrals
setupMM(transport_database, thermo, mode, log_level, trParam);
// do model-specific initialization
tran->initGas(trParam);
}
void TransportFactory::initLiquidTransport(Transport* tran,
thermo_t* thermo,
int log_level)
@ -622,172 +380,6 @@ void TransportFactory::initSolidTransport(Transport* tran,
tran->initSolid(trParam);
}
void TransportFactory::fitCollisionIntegrals(GasTransportParams& tr,
MMCollisionInt& integrals)
{
vector_fp::iterator dptr;
doublereal dstar;
size_t nsp = tr.nsp_;
int mode = tr.mode_;
size_t i, j;
// Chemkin fits to sixth order polynomials
int degree = (mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("tstar_fits\n"
"fits to A*, B*, and C* vs. log(T*).\n"
"These are done only for the required dstar(j,k) values.\n\n");
if (tr.log_level < 3) {
writelog("*** polynomial coefficients not printed (log_level < 3) ***\n");
}
}
for (i = 0; i < nsp; i++) {
for (j = i; j < nsp; j++) {
// Chemkin fits only delta* = 0
if (mode != CK_Mode) {
dstar = tr.delta(i,j);
} else {
dstar = 0.0;
}
// if a fit has already been generated for
// delta* = tr.delta(i,j), then use it. Otherwise,
// make a new fit, and add tr.delta(i,j) to the list
// of delta* values for which fits have been done.
// 'find' returns a pointer to end() if not found
dptr = find(tr.fitlist.begin(), tr.fitlist.end(), dstar);
if (dptr == tr.fitlist.end()) {
vector_fp ca(degree+1), cb(degree+1), cc(degree+1);
vector_fp co22(degree+1);
integrals.fit(degree, dstar,
DATA_PTR(ca), DATA_PTR(cb), DATA_PTR(cc));
integrals.fit_omega22(degree, dstar,
DATA_PTR(co22));
tr.omega22_poly.push_back(co22);
tr.astar_poly.push_back(ca);
tr.bstar_poly.push_back(cb);
tr.cstar_poly.push_back(cc);
tr.poly[i][j] = static_cast<int>(tr.astar_poly.size()) - 1;
tr.fitlist.push_back(dstar);
}
// delta* found in fitlist, so just point to this
// polynomial
else {
tr.poly[i][j] = static_cast<int>((dptr - tr.fitlist.begin()));
}
tr.poly[j][i] = tr.poly[i][j];
}
}
}
void TransportFactory::getTransportData(const ThermoPhase& thermo, const std::vector<const XML_Node*> &xspecies,
XML_Node& log, const std::vector<std::string> &names, GasTransportParams& tr)
{
std::map<std::string, size_t> speciesIndices;
for (size_t i = 0; i < names.size(); i++) {
speciesIndices[names[i]] = i;
}
for (size_t i = 0; i < xspecies.size(); i++) {
const XML_Node& sp = *xspecies[i];
// Find the index for this species in 'names'
size_t j = getValue(speciesIndices, sp["name"], npos);
if (j == npos) {
// Don't need transport data for this species
continue;
}
XML_Node& node = sp.child("transport");
// parameters are converted to SI units before storing
double nAtoms = 0;
size_t kSpec = thermo.speciesIndex(sp["name"]);
for (size_t m = 0; m < thermo.nElements(); m++) {
nAtoms += thermo.nAtoms(kSpec, m);
}
// Molecular geometry; rotational heat capacity / R
XML_Node* geomNode = ctml::getByTitle(node, "geometry");
std::string geom = (geomNode) ? geomNode->value() : "";
if (geom == "atom") {
if (nAtoms != 1) {
throw TransportDBError(i, "invalid geometry. 'atom' specified,"
" but species contains multiple atoms.");
}
tr.crot[j] = 0.0;
} else if (geom == "linear") {
if (nAtoms == 1) {
throw TransportDBError(i, "invalid geometry. 'linear' specified,"
" but species only contains one atom.");
}
tr.crot[j] = 1.0;
} else if (geom == "nonlinear") {
if (nAtoms < 3) {
throw TransportDBError(i, "invalid geometry. 'nonlinear' specified,"
" but species only contains " + fp2str(nAtoms) + " atoms.");
}
tr.crot[j] = 1.5;
} else {
throw TransportDBError(i, "invalid geometry");
}
// Pitzer's acentric factor:
double acentric;
ctml::getOptionalFloat(node, "acentric_factor", acentric);
if (acentric) {
tr.w_ac[j] = acentric;
} /*else {
throw TransportDBError(i, "acentric factor not defined");
}*/
// Well-depth parameter in Kelvin (converted to Joules)
double welldepth = ctml::getFloat(node, "LJ_welldepth");
if (welldepth >= 0.0) {
tr.eps[j] = Boltzmann * welldepth;
} else {
throw TransportDBError(i, "negative well depth");
}
// Lennard-Jones diameter of the molecule, given in Angstroms.
double diam = ctml::getFloat(node, "LJ_diameter");
if (diam > 0.0) {
tr.sigma[j] = 1.e-10 * diam; // A -> m
} else {
throw TransportDBError(i, "negative or zero diameter");
}
// Dipole moment of the molecule.
// Given in Debye (a Debye is 1e-18 statC-m or 3.3356e-30 C-m)
double dipole = ctml::getFloat(node, "dipoleMoment");
if (dipole >= 0.0) {
tr.dipole(j,j) = 1e-21 / lightSpeed * dipole;
tr.polar[j] = (dipole > 0.0);
} else {
throw TransportDBError(i, "negative dipole moment");
}
// Polarizability of the molecule, given in cubic Angstroms.
double polar = ctml::getFloat(node, "polarizability");
if (polar >= 0.0) {
tr.alpha[j] = 1.e-30 * polar; // A^3 -> m^3
} else {
throw TransportDBError(i, "negative polarizability");
}
// Rotational relaxation number. (Number of collisions it takes to
// equilibrate the rotational dofs with the temperature)
double rot = ctml::getFloat(node, "rotRelax");
if (rot >= 0.0) {
tr.zrot[j] = std::max(1.0, rot);
} else {
throw TransportDBError(i, "negative rotation relaxation number");
}
}
}
void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML_Node*> &xspecies,
XML_Node& log,
const std::vector<std::string> &names,
@ -1104,263 +696,6 @@ void TransportFactory::getSolidTransportData(const XML_Node& transportNode,
return;
}
void TransportFactory::fitProperties(GasTransportParams& tr,
MMCollisionInt& integrals)
{
doublereal tstar;
int ndeg = 0;
// number of points to use in generating fit data
const size_t np = 50;
int mode = tr.mode_;
int degree = (mode == CK_Mode ? 3 : 4);
doublereal t, om22;
doublereal dt = (tr.tmax - tr.tmin)/(np-1);
vector_fp tlog(np), spvisc(np), spcond(np);
doublereal val, fit;
vector_fp w(np), w2(np);
// generate array of log(t) values
for (size_t n = 0; n < np; n++) {
t = tr.tmin + dt*n;
tlog[n] = log(t);
}
// vector of polynomial coefficients
vector_fp c(degree + 1), c2(degree + 1);
// fit the pure-species viscosity and thermal conductivity for
// each species
if (DEBUG_MODE_ENABLED && tr.log_level < 2 && m_verbose) {
writelog("*** polynomial coefficients not printed (log_level < 2) ***\n");
}
doublereal sqrt_T, visc, err, relerr,
mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0;
if (DEBUG_MODE_ENABLED && m_verbose) {
writelog("Polynomial fits for viscosity:\n");
if (mode == CK_Mode) {
writelog("log(viscosity) fit to cubic polynomial in log(T)\n");
} else {
writelogf("viscosity/sqrt(T) fit to polynomial of degree "
"%d in log(T)", degree);
}
}
doublereal cp_R, cond, w_RT, f_int, A_factor, B_factor,
c1, cv_rot, cv_int, f_rot, f_trans, om11;
doublereal diffcoeff;
for (size_t k = 0; k < tr.nsp_; k++) {
for (size_t n = 0; n < np; n++) {
t = tr.tmin + dt*n;
tr.thermo->setTemperature(t);
vector_fp cp_R_all(tr.thermo->nSpecies());
tr.thermo->getCp_R_ref(&cp_R_all[0]);
cp_R = cp_R_all[k];
tstar = Boltzmann * t/ tr.eps[k];
sqrt_T = sqrt(t);
om22 = integrals.omega22(tstar, tr.delta(k,k));
om11 = integrals.omega11(tstar, tr.delta(k,k));
// self-diffusion coefficient, without polar
// corrections
diffcoeff = ThreeSixteenths *
sqrt(2.0 * Pi/tr.reducedMass(k,k)) *
pow((Boltzmann * t), 1.5)/
(Pi * tr.sigma[k] * tr.sigma[k] * om11);
// viscosity
visc = FiveSixteenths
* sqrt(Pi * tr.mw[k] * Boltzmann * t / Avogadro) /
(om22 * Pi * tr.sigma[k]*tr.sigma[k]);
// thermal conductivity
w_RT = tr.mw[k]/(GasConstant * t);
f_int = w_RT * diffcoeff/visc;
cv_rot = tr.crot[k];
A_factor = 2.5 - f_int;
B_factor = tr.zrot[k] + TwoOverPi
*(FiveThirds * cv_rot + f_int);
c1 = TwoOverPi * A_factor/B_factor;
cv_int = cp_R - 2.5 - cv_rot;
f_rot = f_int * (1.0 + c1);
f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5);
cond = (visc/tr.mw[k])*GasConstant*(f_trans * 1.5
+ f_rot * cv_rot + f_int * cv_int);
if (mode == CK_Mode) {
spvisc[n] = log(visc);
spcond[n] = log(cond);
w[n] = -1.0;
w2[n] = -1.0;
} else {
// the viscosity should be proportional
// approximately to sqrt(T); therefore,
// visc/sqrt(T) should have only a weak
// temperature dependence. And since the mixture
// rule requires the square root of the
// pure-species viscosity, fit the square root of
// (visc/sqrt(T)) to avoid having to compute
// square roots in the mixture rule.
spvisc[n] = sqrt(visc/sqrt_T);
// the pure-species conductivity scales
// approximately with sqrt(T). Unlike the
// viscosity, there is no reason here to fit the
// square root, since a different mixture rule is
// used.
spcond[n] = cond/sqrt_T;
w[n] = 1.0/(spvisc[n]*spvisc[n]);
w2[n] = 1.0/(spcond[n]*spcond[n]);
}
}
polyfit(np, DATA_PTR(tlog), DATA_PTR(spvisc),
DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c));
polyfit(np, DATA_PTR(tlog), DATA_PTR(spcond),
DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c2));
// evaluate max fit errors for viscosity
for (size_t n = 0; n < np; n++) {
if (mode == CK_Mode) {
val = exp(spvisc[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c)));
} else {
sqrt_T = exp(0.5*tlog[n]);
val = sqrt_T * pow(spvisc[n],2);
fit = sqrt_T * pow(poly4(tlog[n], DATA_PTR(c)),2);
}
err = fit - val;
relerr = err/val;
mxerr = std::max(mxerr, fabs(err));
mxrelerr = std::max(mxrelerr, fabs(relerr));
}
// evaluate max fit errors for conductivity
for (size_t n = 0; n < np; n++) {
if (mode == CK_Mode) {
val = exp(spcond[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c2)));
} else {
sqrt_T = exp(0.5*tlog[n]);
val = sqrt_T * spcond[n];
fit = sqrt_T * poly4(tlog[n], DATA_PTR(c2));
}
err = fit - val;
relerr = err/val;
mxerr_cond = std::max(mxerr_cond, fabs(err));
mxrelerr_cond = std::max(mxrelerr_cond, fabs(relerr));
}
tr.visccoeffs.push_back(c);
tr.condcoeffs.push_back(c2);
if (DEBUG_MODE_ENABLED && tr.log_level >= 2 && m_verbose) {
writelog(tr.thermo->speciesName(k) + ": [" + vec2str(c) + "]\n");
}
}
if (DEBUG_MODE_ENABLED && m_verbose) {
writelogf("Maximum viscosity absolute error: %12.6g\n", mxerr);
writelogf("Maximum viscosity relative error: %12.6g\n", mxrelerr);
writelog("\nPolynomial fits for conductivity:\n");
if (mode == CK_Mode)
writelog("log(conductivity) fit to cubic polynomial in log(T)");
else {
writelogf("conductivity/sqrt(T) fit to "
"polynomial of degree %d in log(T)", degree);
}
if (tr.log_level >= 2)
for (size_t k = 0; k < tr.nsp_; k++) {
writelog(tr.thermo->speciesName(k) + ": [" +
vec2str(tr.condcoeffs[k]) + "]\n");
}
writelogf("Maximum conductivity absolute error: %12.6g\n", mxerr_cond);
writelogf("Maximum conductivity relative error: %12.6g\n", mxrelerr_cond);
// fit the binary diffusion coefficients for each species pair
writelogf("\nbinary diffusion coefficients:\n");
if (mode == CK_Mode)
writelog("log(D) fit to cubic polynomial in log(T)");
else {
writelogf("D/T**(3/2) fit to polynomial of degree %d in log(T)",degree);
}
}
mxerr = 0.0, mxrelerr = 0.0;
vector_fp diff(np + 1);
doublereal eps, sigma;
for (size_t k = 0; k < tr.nsp_; k++) {
for (size_t j = k; j < tr.nsp_; j++) {
for (size_t n = 0; n < np; n++) {
t = tr.tmin + dt*n;
eps = tr.epsilon(j,k);
tstar = Boltzmann * t/eps;
sigma = tr.diam(j,k);
om11 = integrals.omega11(tstar, tr.delta(j,k));
diffcoeff = ThreeSixteenths *
sqrt(2.0 * Pi/tr.reducedMass(k,j)) *
pow((Boltzmann * t), 1.5)/
(Pi * sigma * sigma * om11);
// 2nd order correction
// NOTE: THIS CORRECTION IS NOT APPLIED
doublereal fkj, fjk;
getBinDiffCorrection(t, tr, integrals, k, j, 1.0, 1.0, fkj, fjk);
if (mode == CK_Mode) {
diff[n] = log(diffcoeff);
w[n] = -1.0;
} else {
diff[n] = diffcoeff/pow(t, 1.5);
w[n] = 1.0/(diff[n]*diff[n]);
}
}
polyfit(np, DATA_PTR(tlog), DATA_PTR(diff),
DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c));
doublereal pre;
for (size_t n = 0; n < np; n++) {
if (mode == CK_Mode) {
val = exp(diff[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c)));
} else {
t = exp(tlog[n]);
pre = pow(t, 1.5);
val = pre * diff[n];
fit = pre * poly4(tlog[n], DATA_PTR(c));
}
err = fit - val;
relerr = err/val;
mxerr = std::max(mxerr, fabs(err));
mxrelerr = std::max(mxrelerr, fabs(relerr));
}
tr.diffcoeffs.push_back(c);
if (DEBUG_MODE_ENABLED && tr.log_level >= 2 && m_verbose) {
writelog(tr.thermo->speciesName(k) + "__" +
tr.thermo->speciesName(j) + ": [" + vec2str(c) + "]\n");
}
}
}
if (DEBUG_MODE_ENABLED && m_verbose) {
writelogf("Maximum binary diffusion coefficient absolute error:"
" %12.6g\n", mxerr);
writelogf("Maximum binary diffusion coefficient relative error:"
"%12.6g", mxrelerr);
}
}
Transport* newTransportMgr(const std::string& transportModel, thermo_t* thermo, int loglevel, TransportFactory* f, int ndim)
{