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

[Transport] Construct gas phase Transport without GasTransportParams

Browse source 
Since most of the members of class GasTransportParams were copied to
MixTransport or MultiTransport, they have instead been made members of
GasTransport and populated directly, making class GasTransportParams
unnecessary.
This commit is contained in:
Ray Speth 2015-01-13 23:07:00 +00:00
parent 44a99de4e0
commit fdf7db716e
13 changed files with 396 additions and 412 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -107,23 +107,11 @@ 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);
virtual bool initGas(GasTransportParams& tr);
virtual void update_T();
virtual void update_C() = 0;
@ -159,35 +147,17 @@ protected:
//! 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 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
* Uses polynomial fits to Monchick & Mason collision integrals.
*/
void setupMM(thermo_t* thermo, int mode, int log_level,
GasTransportParams& tr);
void setupMM();
//! 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 tr Output object containing the transport parameters for
* the species listed in names (in the order of their
* listing in names).
* species names.
*/
void getTransportData(const ThermoPhase& thermo, GasTransportParams& tr);
void getTransportData();
//! Corrections for polar-nonpolar binary diffusion coefficients
/*!
@ -198,22 +168,17 @@ protected:
*
* @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,
void makePolarCorrections(size_t i, size_t j, 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);
void fitCollisionIntegrals(MMCollisionInt& integrals);
//! Generate polynomial fits to the viscosity, conductivity, and
//! the binary diffusion coefficients
@ -235,11 +200,9 @@ protected:
* 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);
void fitProperties(MMCollisionInt& integrals);
//! Second-order correction to the binary diffusion coefficients
/*!
@ -253,7 +216,6 @@ protected:
* 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
@ -264,8 +226,7 @@ protected:
*
* @note This method is not used currently.
*/
void getBinDiffCorrection(doublereal t, const GasTransportParams& tr,
MMCollisionInt& integrals, size_t k,
void getBinDiffCorrection(doublereal t, MMCollisionInt& integrals, size_t k,
size_t j, doublereal xk, doublereal xj,
doublereal& fkj, doublereal& fjk);
@ -360,9 +321,9 @@ protected:
//! Polynomial fits to the binary diffusivity of each species
/*!
* m_diffcoeff[ic] is vector of polynomial coefficients for species i species j
* that fits the binary diffusion coefficient. The relationship between i
* j and ic is determined from the following algorithm:
* m_diffcoeff[ic] is vector of polynomial coefficients for species i
* species j that fits the binary diffusion coefficient. The relationship
* between i j and ic is determined from the following algorithm:
*
* int ic = 0;
* for (i = 0; i < m_nsp; i++) {
@ -377,9 +338,149 @@ protected:
//! the current temperature Size is nsp x nsp.
DenseMatrix m_bdiff;
//! Boolean indicating whether to turn on verbose printing during
//! initialization
bool m_verbose;
//! temperature fits of the heat conduction
/*!
* Dimensions are number of species (nsp) polynomial order of the collision
* integral fit (degree+1).
*/
std::vector<vector_fp> m_condcoeffs;
//! Indices for the (i,j) interaction in collision integral fits
/*!
* m_poly[i][j] contains the index for (i,j) interactions in
* #m_omega22_poly, #m_astar_poly, #m_bstar_poly, and #m_cstar_poly.
*/
std::vector<vector_int> m_poly;
//! Fit for omega22 collision integral
/*!
* m_omega22_poly[m_poly[i][j]] is the vector of polynomial coefficients
* (length degree+1) for the collision integral fit for the species pair
* (i,j).
*/
std::vector<vector_fp> m_omega22_poly;
//! Fit for astar collision integral
/*!
* m_astar_poly[m_poly[i][j]] is the vector of polynomial coefficients
* (length degree+1) for the collision integral fit for the species pair
* (i,j).
*/
std::vector<vector_fp> m_astar_poly;
//! Fit for bstar collision integral
/*!
* m_bstar_poly[m_poly[i][j]] is the vector of polynomial coefficients
* (length degree+1) for the collision integral fit for the species pair
* (i,j).
*/
std::vector<vector_fp> m_bstar_poly;
//! Fit for cstar collision integral
/*!
* m_bstar_poly[m_poly[i][j]] is the vector of polynomial coefficients
* (length degree+1) for the collision integral fit for the species pair
* (i,j).
*/
std::vector<vector_fp> m_cstar_poly;
//! Rotational relaxation number for each species
/*!
* length is the number of species in the phase. units are dimensionless
*/
vector_fp m_zrot;
//! Dimensionless rotational heat capacity of each species
/*!
* These values are 0, 1 and 1.5 for single-molecule, linear, and nonlinear
* species respectively length is the number of species in the phase.
* Dimensionless (Cr / R)
*/
vector_fp m_crot;
//! Vector of booleans indicating whether a species is a polar molecule
/*!
* Length is nsp
*/
std::vector<bool> m_polar;
//! Polarizability of each species in the phase
/*!
* Length = nsp. Units = m^3
*/
vector_fp m_alpha;
//! Lennard-Jones well-depth of the species in the current phase
/*!
* length is the number of species in the phase. Units are Joules (Note this
* is not Joules/kmol) (note, no kmol -> this is a per molecule amount)
*/
vector_fp m_eps;
//! Lennard-Jones diameter of the species in the current phase
/*!
* length is the number of species in the phase
* units are in meters.
*/
vector_fp m_sigma;
//! This is the reduced mass of the interaction between species i and j
/*!
* reducedMass(i,j) = mw[i] * mw[j] / (Avogadro * (mw[i] + mw[j]));
*
* Units are kg (note, no kmol -> this is a per molecule amount)
*
* Length nsp * nsp. This is a symmetric matrix
*/
DenseMatrix m_reducedMass;
//! hard-sphere diameter for (i,j) collision
/*!
* diam(i,j) = 0.5*(sigma[i] + sigma[j]);
* Units are m (note, no kmol -> this is a per molecule amount)
*
* Length nsp * nsp. This is a symmetric matrix.
*/
DenseMatrix m_diam;
//! The effective well depth for (i,j) collisions
/*!
* epsilon(i,j) = sqrt(eps[i]*eps[j]);
* Units are Joules (note, no kmol -> this is a per molecule amount)
*
* Length nsp * nsp. This is a symmetric matrix.
*/
DenseMatrix m_epsilon;
//! The effective dipole moment for (i,j) collisions
/*!
* Given `dipoleMoment` in Debye (a Debye is 3.335e-30 C-m):
*
* dipole(i,i) = 1.e-21 / lightSpeed * dipoleMoment;
* dipole(i,j) = sqrt(dipole(i,i) * dipole(j,j));
* (note, no kmol -> this is a per molecule amount)
*
* Length nsp * nsp. This is a symmetric matrix.
*/
DenseMatrix m_dipole;
//! Reduced dipole moment of the interaction between two species
/*!
* This is the reduced dipole moment of the interaction between two species
* 0.5 * dipole(i,j)^2 / (4 * Pi * epsilon_0 * epsilon(i,j) * d^3);
*
* Length nsp * nsp .This is a symmetric matrix
*/
DenseMatrix m_delta;
//! Pitzer acentric factor
/*!
* Length is the number of species in the phase. Dimensionless.
*/
vector_fp m_w_ac;
//! Level of verbose printing during initialization
int m_log_level;
};
} // namespace Cantera

8
include/cantera/transport/HighPressureGasTransport.h
View file

@ -14,8 +14,6 @@
namespace Cantera
{
class GasTransportParams;
//! Class MultiTransport implements transport properties for
//! high pressure gas mixtures.
/*!
@ -68,12 +66,6 @@ public:
virtual doublereal viscosity();
//! Initialize the transport operator with parameters from GasTransportParams object
/*!
* @param tr input GasTransportParams object
*/
virtual bool initGas(GasTransportParams& tr);
friend class TransportFactory;
protected:

28
include/cantera/transport/MixTransport.h
View file

@ -157,17 +157,7 @@ public:
size_t ldx, const doublereal* const grad_X,
size_t ldf, doublereal* const fluxes);
//! Initialize the transport object
/*!
* Here we change all of the internal dimensions to be sufficient.
* We get the object ready to do property evaluations.
*
* @param tr Transport parameters for all of the species
* in the phase.
*/
virtual bool initGas(GasTransportParams& tr);
friend class TransportFactory;
virtual void init(thermo_t* thermo, int mode=0, int log_level=0);
private:
@ -185,13 +175,6 @@ private:
void updateCond_T();
private:
//! Polynomial fits to the thermal conductivity of each species
/*!
* m_condcoeffs[k] is vector of polynomial coefficients for species k
* that fits the thermal conductivity
*/
std::vector<vector_fp> m_condcoeffs;
//! vector of species thermal conductivities (W/m /K)
/*!
* These are used in wilke's rule to calculate the viscosity of the
@ -210,14 +193,7 @@ private:
//! Update boolean for the mixture rule for the mixture thermal conductivity
bool m_condmix_ok;
public:
vector_fp m_eps;
vector_fp m_sigma;
vector_fp m_alpha;
DenseMatrix m_dipole;
vector_fp m_zrot;
vector_fp m_crot;
private:
//! Debug flag - turns on more printing
bool m_debug;
};

24
include/cantera/transport/MultiTransport.h
View file

@ -115,13 +115,7 @@ public:
const doublereal* state2, doublereal delta,
doublereal* fluxes);
//! Initialize the transport operator with parameters from GasTransportParams object
/*!
* @param tr input GasTransportParams object
*/
virtual bool initGas(GasTransportParams& tr);
friend class TransportFactory;
virtual void init(ThermoPhase* thermo, int mode=0, int log_level=0);
protected:
//! Update basic temperature-dependent quantities if the temperature has changed.
@ -136,13 +130,6 @@ protected:
doublereal m_thermal_tlast;
// property values
std::vector<std::vector<int> > m_poly;
std::vector<vector_fp> m_astar_poly;
std::vector<vector_fp> m_bstar_poly;
std::vector<vector_fp> m_cstar_poly;
std::vector<vector_fp> m_om22_poly;
//! Dense matrix for astar
DenseMatrix m_astar;
@ -155,17 +142,8 @@ protected:
//! Dense matrix for omega22
DenseMatrix m_om22;
public:
vector_fp m_crot;
vector_fp m_cinternal;
vector_fp m_zrot;
vector_fp m_eps;
vector_fp m_sigma;
vector_fp m_alpha;
vector_fp m_w_ac;
DenseMatrix m_dipole;
protected:
vector_fp m_sqrt_eps_k;
DenseMatrix m_log_eps_k;
vector_fp m_frot_298;

2
include/cantera/transport/PecosTransport.h
View file

@ -13,6 +13,8 @@
namespace Cantera
{
class GasTransportParams;
/**
* Class PecosTransport implements mixture-averaged transport
* properties for ideal gas mixtures.

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

@ -26,7 +26,6 @@
namespace Cantera
{
class GasTransportParams;
class LiquidTransportParams;
class SolidTransportData;
@ -713,26 +712,24 @@ public:
return m_velocityBasis;
}
friend class TransportFactory;
protected:
/**
* @name Transport manager construction
* These methods are used internally during construction.
* These methods are used during construction.
* @{
*/
//! Called by TransportFactory to set parameters.
//! Initialize a transport manager
/*!
* This is called by classes that use the gas phase parameter
* list to initialize themselves.
* This routine sets up a transport manager. It calculates the collision
* integrals and populates species-dependent data structures.
*
* @param tr Reference to the parameter list that will be used
* to initialize the class
* @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 bool initGas(GasTransportParams& tr) {
throw NotImplementedError("Transport::initGas");
virtual void init(thermo_t* thermo, int mode=0, int log_level=0) {
throw NotImplementedError("Transport::init");
}
//! Called by TransportFactory to set parameters.
@ -747,7 +744,6 @@ protected:
throw NotImplementedError("Transport::initLiquid");
}
public:
//! Called by TransportFactory to set parameters.
/*!
* This is called by classes that use the solid phase parameter

2
include/cantera/transport/TransportParams.h
View file

@ -58,6 +58,8 @@ public:
//! gases with a kinetic theory of gases derived transport model.
/*!
* This structure is used by TransportFactory object.
* @deprecated Unused. Mostly merged into class GasTransport. This class will be
* removed after Cantera 2.2.
*/
class GasTransportParams : public TransportParams
{

436
src/transport/GasTransport.cpp
View file

@ -1,6 +1,5 @@
//! @file GasTransport.cpp
#include "cantera/transport/GasTransport.h"
#include "cantera/transport/TransportParams.h"
#include "MMCollisionInt.h"
#include "cantera/base/stringUtils.h"
#include "cantera/numerics/polyfit.h"
@ -40,7 +39,7 @@ GasTransport::GasTransport(ThermoPhase* thermo) :
m_t32(0.0),
m_diffcoeffs(0),
m_bdiff(0, 0),
m_verbose(false)
m_log_level(0)
{
}
@ -70,7 +69,7 @@ GasTransport::GasTransport(const GasTransport& right) :
m_t32(0.0),
m_diffcoeffs(0),
m_bdiff(0, 0),
m_verbose(false)
m_log_level(0)
{
}
@ -100,54 +99,28 @@ 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;
m_condcoeffs = right.m_condcoeffs;
m_poly = right.m_poly;
m_omega22_poly = right.m_omega22_poly;
m_astar_poly = right.m_astar_poly;
m_bstar_poly = right.m_bstar_poly;
m_cstar_poly = right.m_cstar_poly;
m_zrot = right.m_zrot;
m_polar = right.m_polar;
m_alpha = right.m_alpha;
m_eps = right.m_eps;
m_sigma = right.m_sigma;
m_reducedMass = right.m_reducedMass;
m_diam = right.m_diam;
m_epsilon = right.m_epsilon;
m_dipole = right.m_dipole;
m_delta = right.m_delta;
m_w_ac = right.m_w_ac;
m_log_level = right.m_log_level;
return *this;
}
bool GasTransport::initGas(GasTransportParams& tr)
{
// constant mixture attributes
m_thermo = tr.thermo;
m_nsp = m_thermo->nSpecies();
// copy polynomials and parameters into local storage
m_visccoeffs = tr.visccoeffs;
m_diffcoeffs = tr.diffcoeffs;
m_mode = tr.mode_;
m_molefracs.resize(m_nsp);
m_spwork.resize(m_nsp);
m_visc.resize(m_nsp);
m_phi.resize(m_nsp, m_nsp, 0.0);
m_bdiff.resize(m_nsp, m_nsp);
// 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());
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
for (size_t j = 0; j < m_nsp; j++) {
for (size_t 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_sqvisc.resize(m_nsp);
// set flags all false
m_visc_ok = false;
m_viscwt_ok = false;
m_spvisc_ok = false;
m_bindiff_ok = false;
return true;
}
void GasTransport::update_T(void)
{
double T = m_thermo->temperature();
@ -385,150 +358,160 @@ 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;
}
m_thermo = thermo;
m_nsp = m_thermo->nSpecies();
m_mode = mode;
m_log_level = log_level;
// set up Monchick and Mason collision integrals
setupMM(thermo, mode, log_level, trParam);
// do model-specific initialization
initGas(trParam);
setupMM();
m_molefracs.resize(m_nsp);
m_spwork.resize(m_nsp);
m_visc.resize(m_nsp);
m_sqvisc.resize(m_nsp);
m_phi.resize(m_nsp, m_nsp, 0.0);
m_bdiff.resize(m_nsp, m_nsp);
// make a local copy of the molecular weights
m_mw.assign(m_thermo->molecularWeights().begin(),
m_thermo->molecularWeights().end());
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
for (size_t j = 0; j < m_nsp; j++) {
for (size_t 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]);
}
}
// set flags all false
m_visc_ok = false;
m_viscwt_ok = false;
m_spvisc_ok = false;
m_bindiff_ok = false;
}
void GasTransport::setupMM(thermo_t* thermo, int mode, int log_level,
GasTransportParams& tr)
void GasTransport::setupMM()
{
// constant mixture attributes
tr.thermo = thermo;
tr.nsp_ = tr.thermo->nSpecies();
size_t nsp = tr.nsp_;
m_epsilon.resize(m_nsp, m_nsp, 0.0);
m_delta.resize(m_nsp, m_nsp, 0.0);
m_reducedMass.resize(m_nsp, m_nsp, 0.0);
m_dipole.resize(m_nsp, m_nsp, 0.0);
m_diam.resize(m_nsp, m_nsp, 0.0);
m_crot.resize(m_nsp);
m_zrot.resize(m_nsp);
m_polar.resize(m_nsp, false);
m_alpha.resize(m_nsp, 0.0);
m_poly.resize(m_nsp);
m_sigma.resize(m_nsp);
m_eps.resize(m_nsp);
m_w_ac.resize(m_nsp);
tr.tmin = thermo->minTemp();
tr.tmax = thermo->maxTemp();
tr.mw.resize(nsp);
tr.log_level = log_level;
const vector_fp& mw = m_thermo->molecularWeights();
getTransportData();
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);
getTransportData(*thermo, tr);
for (size_t i = 0; i < nsp; i++) {
tr.poly[i].resize(nsp);
for (size_t i = 0; i < m_nsp; i++) {
m_poly[i].resize(m_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++) {
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = i; j < m_nsp; j++) {
// the reduced mass
tr.reducedMass(i,j) = tr.mw[i] * tr.mw[j] / (Avogadro * (tr.mw[i] + tr.mw[j]));
m_reducedMass(i,j) = mw[i] * mw[j] / (Avogadro * (mw[i] + mw[j]));
// hard-sphere diameter for (i,j) collisions
tr.diam(i,j) = 0.5*(tr.sigma[i] + tr.sigma[j]);
m_diam(i,j) = 0.5*(m_sigma[i] + m_sigma[j]);
// the effective well depth for (i,j) collisions
tr.epsilon(i,j) = sqrt(tr.eps[i]*tr.eps[j]);
m_epsilon(i,j) = sqrt(m_eps[i]*m_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));
tstar_min = std::min(tstar_min, Boltzmann * m_thermo->minTemp()/m_epsilon(i,j));
tstar_max = std::max(tstar_max, Boltzmann * m_thermo->maxTemp()/m_epsilon(i,j));
// the effective dipole moment for (i,j) collisions
tr.dipole(i,j) = sqrt(tr.dipole(i,i)*tr.dipole(j,j));
m_dipole(i,j) = sqrt(m_dipole(i,i)*m_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);
double d = m_diam(i,j);
m_delta(i,j) = 0.5 * m_dipole(i,j)*m_dipole(i,j)
/ (4 * Pi * epsilon_0 * m_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;
makePolarCorrections(i, j, f_eps, f_sigma);
m_diam(i,j) *= f_sigma;
m_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);
m_reducedMass(j,i) = m_reducedMass(i,j);
m_diam(j,i) = m_diam(i,j);
m_epsilon(j,i) = m_epsilon(i,j);
m_dipole(j,i) = m_dipole(i,j);
m_delta(j,i) = m_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) {
if (m_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) {
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** collision_integrals ***\n");
}
MMCollisionInt integrals;
integrals.init(tstar_min, tstar_max, log_level);
fitCollisionIntegrals(tr, integrals);
if (DEBUG_MODE_ENABLED && m_verbose) {
integrals.init(tstar_min, tstar_max, m_log_level);
fitCollisionIntegrals(integrals);
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** end of collision_integrals ***\n");
}
// make polynomial fits
if (DEBUG_MODE_ENABLED && m_verbose) {
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** property fits ***\n");
}
fitProperties(tr, integrals);
if (DEBUG_MODE_ENABLED && m_verbose) {
fitProperties(integrals);
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** end of property fits ***\n");
}
}
void GasTransport::getTransportData(const ThermoPhase& thermo,
GasTransportParams& tr)
void GasTransport::getTransportData()
{
for (size_t k = 0; k < thermo.nSpecies(); k++) {
const Species& s = thermo.species(thermo.speciesName(k));
for (size_t k = 0; k < m_thermo->nSpecies(); k++) {
const Species& s = m_thermo->species(m_thermo->speciesName(k));
const GasTransportData& sptran =
dynamic_cast<GasTransportData&>(*s.transport.get());
if (sptran.geometry == "atom") {
tr.crot[k] = 0.0;
m_crot[k] = 0.0;
} else if (sptran.geometry == "linear") {
tr.crot[k] = 1.0;
m_crot[k] = 1.0;
} else if (sptran.geometry == "nonlinear") {
tr.crot[k] = 1.5;
m_crot[k] = 1.5;
}
tr.sigma[k] = sptran.diameter;
tr.eps[k] = sptran.well_depth;
tr.dipole(k,k) = sptran.dipole;
tr.polar[k] = (sptran.dipole > 0);
tr.alpha[k] = sptran.polarizability;
tr.zrot[k] = sptran.rotational_relaxation;
tr.w_ac[k] = sptran.acentric_factor;
m_sigma[k] = sptran.diameter;
m_eps[k] = sptran.well_depth;
m_dipole(k,k) = sptran.dipole;
m_polar[k] = (sptran.dipole > 0);
m_alpha[k] = sptran.polarizability;
m_zrot[k] = sptran.rotational_relaxation;
m_w_ac[k] = sptran.acentric_factor;
}
}
void GasTransport::makePolarCorrections(size_t i, size_t j,
const GasTransportParams& tr, doublereal& f_eps, doublereal& f_sigma)
doublereal& f_eps, doublereal& f_sigma)
{
// no correction if both are nonpolar, or both are polar
if (tr.polar[i] == tr.polar[j]) {
if (m_polar[i] == m_polar[j]) {
f_eps = 1.0;
f_sigma = 1.0;
return;
@ -537,94 +520,89 @@ void GasTransport::makePolarCorrections(size_t i, size_t j,
// 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 kp = (m_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]);
d3np = pow(m_sigma[knp],3);
d3p = pow(m_sigma[kp],3);
alpha_star = m_alpha[knp]/d3np;
mu_p_star = m_dipole(kp,kp)/sqrt(4 * Pi * epsilon_0 * d3p * m_eps[kp]);
xi = 1.0 + 0.25 * alpha_star * mu_p_star * mu_p_star *
sqrt(tr.eps[kp]/tr.eps[knp]);
sqrt(m_eps[kp]/m_eps[knp]);
f_sigma = pow(xi, -1.0/6.0);
f_eps = xi*xi;
}
void GasTransport::fitCollisionIntegrals(GasTransportParams& tr,
MMCollisionInt& integrals)
void GasTransport::fitCollisionIntegrals(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) {
int degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
if (DEBUG_MODE_ENABLED && m_log_level) {
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) {
if (m_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++) {
vector_fp fitlist;
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = i; j < m_nsp; j++) {
// Chemkin fits only delta* = 0
if (mode != CK_Mode) {
dstar = tr.delta(i,j);
if (m_mode != CK_Mode) {
dstar = m_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
// if a fit has already been generated for delta* = m_delta(i,j),
// then use it. Otherwise, make a new fit, and add m_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::iterator dptr = find(fitlist.begin(), fitlist.end(), dstar);
if (dptr == 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);
m_omega22_poly.push_back(co22);
m_astar_poly.push_back(ca);
m_bstar_poly.push_back(cb);
m_cstar_poly.push_back(cc);
m_poly[i][j] = static_cast<int>(m_astar_poly.size()) - 1;
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()));
m_poly[i][j] = static_cast<int>((dptr - fitlist.begin()));
}
tr.poly[j][i] = tr.poly[i][j];
m_poly[j][i] = m_poly[i][j];
}
}
}
void GasTransport::fitProperties(GasTransportParams& tr,
MMCollisionInt& integrals)
void GasTransport::fitProperties(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);
int degree = (m_mode == CK_Mode ? 3 : 4);
double dt = (tr.tmax - tr.tmin)/(np-1);
double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(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;
double t = m_thermo->minTemp() + dt*n;
tlog[n] = log(t);
}
@ -632,15 +610,15 @@ void GasTransport::fitProperties(GasTransportParams& tr,
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) {
if (DEBUG_MODE_ENABLED && m_log_level && m_log_level < 2) {
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) {
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("Polynomial fits for viscosity:\n");
if (mode == CK_Mode) {
if (m_mode == CK_Mode) {
writelog("log(viscosity) fit to cubic polynomial in log(T)\n");
} else {
writelogf("viscosity/sqrt(T) fit to polynomial of degree "
@ -651,47 +629,48 @@ void GasTransport::fitProperties(GasTransportParams& tr,
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++) {
const vector_fp& mw = m_thermo->molecularWeights();
for (size_t k = 0; k < m_nsp; k++) {
for (size_t n = 0; n < np; n++) {
double t = tr.tmin + dt*n;
double t = m_thermo->minTemp() + dt*n;
tr.thermo->setTemperature(t);
vector_fp cp_R_all(tr.thermo->nSpecies());
tr.thermo->getCp_R_ref(&cp_R_all[0]);
m_thermo->setTemperature(t);
vector_fp cp_R_all(m_thermo->nSpecies());
m_thermo->getCp_R_ref(&cp_R_all[0]);
cp_R = cp_R_all[k];
double tstar = Boltzmann * t/ tr.eps[k];
double tstar = Boltzmann * t/ m_eps[k];
sqrt_T = sqrt(t);
double om22 = integrals.omega22(tstar, tr.delta(k,k));
om11 = integrals.omega11(tstar, tr.delta(k,k));
double om22 = integrals.omega22(tstar, m_delta(k,k));
om11 = integrals.omega11(tstar, m_delta(k,k));
// self-diffusion coefficient, without polar corrections
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/tr.reducedMass(k,k)) *
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,k)) *
pow((Boltzmann * t), 1.5)/
(Pi * tr.sigma[k] * tr.sigma[k] * om11);
(Pi * m_sigma[k] * m_sigma[k] * om11);
// viscosity
visc = FiveSixteenths
* sqrt(Pi * tr.mw[k] * Boltzmann * t / Avogadro) /
(om22 * Pi * tr.sigma[k]*tr.sigma[k]);
* sqrt(Pi * mw[k] * Boltzmann * t / Avogadro) /
(om22 * Pi * m_sigma[k]*m_sigma[k]);
// thermal conductivity
w_RT = tr.mw[k]/(GasConstant * t);
w_RT = mw[k]/(GasConstant * t);
f_int = w_RT * diffcoeff/visc;
cv_rot = tr.crot[k];
cv_rot = m_crot[k];
A_factor = 2.5 - f_int;
B_factor = tr.zrot[k] + 2.0/Pi * (5.0/3.0 * cv_rot + f_int);
B_factor = m_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
cond = (visc/mw[k])*GasConstant*(f_trans * 1.5
+ f_rot * cv_rot + f_int * cv_int);
if (mode == CK_Mode) {
if (m_mode == CK_Mode) {
spvisc[n] = log(visc);
spcond[n] = log(cond);
w[n] = -1.0;
@ -721,7 +700,7 @@ void GasTransport::fitProperties(GasTransportParams& tr,
// evaluate max fit errors for viscosity
for (size_t n = 0; n < np; n++) {
double val, fit;
if (mode == CK_Mode) {
if (m_mode == CK_Mode) {
val = exp(spvisc[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c)));
} else {
@ -738,7 +717,7 @@ void GasTransport::fitProperties(GasTransportParams& tr,
// evaluate max fit errors for conductivity
for (size_t n = 0; n < np; n++) {
double val, fit;
if (mode == CK_Mode) {
if (m_mode == CK_Mode) {
val = exp(spcond[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c2)));
} else {
@ -751,35 +730,35 @@ void GasTransport::fitProperties(GasTransportParams& tr,
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);
m_visccoeffs.push_back(c);
m_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_log_level >= 2) {
writelog(m_thermo->speciesName(k) + ": [" + vec2str(c) + "]\n");
}
}
if (DEBUG_MODE_ENABLED && m_verbose) {
if (DEBUG_MODE_ENABLED && m_log_level) {
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)
if (m_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");
if (m_log_level >= 2)
for (size_t k = 0; k < m_nsp; k++) {
writelog(m_thermo->speciesName(k) + ": [" +
vec2str(m_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)
if (m_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);
@ -789,25 +768,25 @@ void GasTransport::fitProperties(GasTransportParams& tr,
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 k = 0; k < m_nsp; k++) {
for (size_t j = k; j < m_nsp; j++) {
for (size_t n = 0; n < np; n++) {
double t = tr.tmin + dt*n;
eps = tr.epsilon(j,k);
double t = m_thermo->minTemp() + dt*n;
eps = m_epsilon(j,k);
double tstar = Boltzmann * t/eps;
sigma = tr.diam(j,k);
om11 = integrals.omega11(tstar, tr.delta(j,k));
sigma = m_diam(j,k);
om11 = integrals.omega11(tstar, m_delta(j,k));
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/tr.reducedMass(k,j)) *
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_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);
getBinDiffCorrection(t, integrals, k, j, 1.0, 1.0, fkj, fjk);
if (mode == CK_Mode) {
if (m_mode == CK_Mode) {
diff[n] = log(diffcoeff);
w[n] = -1.0;
} else {
@ -820,7 +799,7 @@ void GasTransport::fitProperties(GasTransportParams& tr,
for (size_t n = 0; n < np; n++) {
double val, fit;
if (mode == CK_Mode) {
if (m_mode == CK_Mode) {
val = exp(diff[n]);
fit = exp(poly3(tlog[n], DATA_PTR(c)));
} else {
@ -834,14 +813,14 @@ void GasTransport::fitProperties(GasTransportParams& tr,
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");
m_diffcoeffs.push_back(c);
if (DEBUG_MODE_ENABLED && m_log_level >= 2) {
writelog(m_thermo->speciesName(k) + "__" +
m_thermo->speciesName(j) + ": [" + vec2str(c) + "]\n");
}
}
}
if (DEBUG_MODE_ENABLED && m_verbose) {
if (DEBUG_MODE_ENABLED && m_log_level) {
writelogf("Maximum binary diffusion coefficient absolute error:"
" %12.6g\n", mxerr);
writelogf("Maximum binary diffusion coefficient relative error:"
@ -849,33 +828,32 @@ void GasTransport::fitProperties(GasTransportParams& tr,
}
}
void GasTransport::getBinDiffCorrection(double t,
const GasTransportParams& tr, MMCollisionInt& integrals,
void GasTransport::getBinDiffCorrection(double t, 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 w1 = m_thermo->molecularWeight(k);
double w2 = m_thermo->molecularWeight(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 sig1 = m_sigma[k];
double sig2 = m_sigma[j];
double sig12 = 0.5*(m_sigma[k] + m_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 tstar1 = Boltzmann * t / m_eps[k];
double tstar2 = Boltzmann * t / m_eps[j];
double tstar12 = Boltzmann * t / sqrt(m_eps[k] * m_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 om22_1 = integrals.omega22(tstar1, m_delta(k,k));
double om22_2 = integrals.omega22(tstar2, m_delta(j,j));
double om11_12 = integrals.omega11(tstar12, m_delta(k,j));
double astar_12 = integrals.astar(tstar12, m_delta(k,j));
double bstar_12 = integrals.bstar(tstar12, m_delta(k,j));
double cstar_12 = integrals.cstar(tstar12, m_delta(k,j));
double cnst = sigratio * sqrt(2.0*w2/wsum) * 2.0 * w1*w1/(wsum * w2);
double p1 = cnst * om22_1 / om11_12;

9
src/transport/HighPressureGasTransport.cpp
View file

@ -32,15 +32,6 @@ HighPressureGasTransport::HighPressureGasTransport(thermo_t* thermo)
{
}
bool HighPressureGasTransport::initGas(GasTransportParams& tr)
{
MultiTransport::initGas(tr);
// copy parameters into local storage
m_w_ac = tr.w_ac;
return true;
}
double HighPressureGasTransport::thermalConductivity()
{
// Method of Ely and Hanley:

20
src/transport/MixTransport.cpp
View file

@ -5,7 +5,6 @@
// copyright 2001 California Institute of Technology
#include "cantera/transport/MixTransport.h"
#include "cantera/transport/TransportParams.h"
#include "cantera/base/stringUtils.h"
using namespace std;
@ -13,7 +12,6 @@ using namespace std;
namespace Cantera
{
MixTransport::MixTransport() :
m_condcoeffs(0),
m_cond(0),
m_lambda(0.0),
m_spcond_ok(false),
@ -24,7 +22,6 @@ MixTransport::MixTransport() :
MixTransport::MixTransport(const MixTransport& right) :
GasTransport(right),
m_condcoeffs(0),
m_cond(0),
m_lambda(0.0),
m_spcond_ok(false),
@ -41,7 +38,6 @@ MixTransport& MixTransport::operator=(const MixTransport& right)
}
GasTransport::operator=(right);
m_condcoeffs = right.m_condcoeffs;
m_cond = right.m_cond;
m_lambda = right.m_lambda;
m_spcond_ok = right.m_spcond_ok;
@ -56,27 +52,15 @@ Transport* MixTransport::duplMyselfAsTransport() const
return new MixTransport(*this);
}
bool MixTransport::initGas(GasTransportParams& tr)
void MixTransport::init(ThermoPhase* thermo, int mode, int log_level)
{
GasTransport::initGas(tr);
m_eps = tr.eps;
m_sigma = tr.sigma;
m_alpha = tr.alpha;
m_dipole = tr.dipole;
m_zrot = tr.zrot;
m_crot = tr.crot;
// copy polynomials and parameters into local storage
m_condcoeffs = tr.condcoeffs;
GasTransport::init(thermo, mode, log_level);
m_cond.resize(m_nsp);
// set flags all false
m_spcond_ok = false;
m_condmix_ok = false;
return true;
}
void MixTransport::getMobilities(doublereal* const mobil)

32
src/transport/MultiTransport.cpp
View file

@ -8,7 +8,6 @@
*/
#include "cantera/transport/MultiTransport.h"
#include "cantera/transport/TransportParams.h"
#include "cantera/thermo/IdealGasPhase.h"
#include "cantera/base/stringUtils.h"
@ -40,23 +39,9 @@ MultiTransport::MultiTransport(thermo_t* thermo)
{
}
bool MultiTransport::initGas(GasTransportParams& tr)
void MultiTransport::init(ThermoPhase* thermo, int mode, int log_level)
{
GasTransport::initGas(tr);
// copy polynomials and parameters into local storage
m_poly = tr.poly;
m_astar_poly = tr.astar_poly;
m_bstar_poly = tr.bstar_poly;
m_cstar_poly = tr.cstar_poly;
m_om22_poly = tr.omega22_poly;
m_zrot = tr.zrot;
m_crot = tr.crot;
m_eps = tr.eps;
m_sigma = tr.sigma;
m_alpha = tr.alpha;
m_dipole = tr.dipole;
m_zrot = tr.zrot;
GasTransport::init(thermo, mode, log_level);
// the L matrix
m_Lmatrix.resize(3*m_nsp, 3*m_nsp);
@ -92,7 +77,7 @@ bool MultiTransport::initGas(GasTransportParams& tr)
// int j;
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = i; j < m_nsp; j++) {
m_log_eps_k(i,j) = log(tr.epsilon(i,j)/Boltzmann);
m_log_eps_k(i,j) = log(m_epsilon(i,j)/Boltzmann);
m_log_eps_k(j,i) = m_log_eps_k(i,j);
}
}
@ -103,12 +88,9 @@ bool MultiTransport::initGas(GasTransportParams& tr)
const doublereal kb298 = Boltzmann * 298.0;
m_sqrt_eps_k.resize(m_nsp);
for (size_t k = 0; k < m_nsp; k++) {
m_sqrt_eps_k[k] = sqrt(tr.eps[k]/Boltzmann);
m_frot_298[k] = Frot(tr.eps[k]/kb298,
m_sqrt_eps_k[k]/sq298);
m_sqrt_eps_k[k] = sqrt(m_eps[k]/Boltzmann);
m_frot_298[k] = Frot(m_eps[k]/kb298, m_sqrt_eps_k[k]/sq298);
}
return true;
}
doublereal MultiTransport::thermalConductivity()
@ -499,12 +481,12 @@ void MultiTransport::updateThermal_T()
z = m_logt - m_log_eps_k(i,j);
ipoly = m_poly[i][j];
if (m_mode == CK_Mode) {
m_om22(i,j) = poly6(z, DATA_PTR(m_om22_poly[ipoly]));
m_om22(i,j) = poly6(z, DATA_PTR(m_omega22_poly[ipoly]));
m_astar(i,j) = poly6(z, DATA_PTR(m_astar_poly[ipoly]));
m_bstar(i,j) = poly6(z, DATA_PTR(m_bstar_poly[ipoly]));
m_cstar(i,j) = poly6(z, DATA_PTR(m_cstar_poly[ipoly]));
} else {
m_om22(i,j) = poly8(z, DATA_PTR(m_om22_poly[ipoly]));
m_om22(i,j) = poly8(z, DATA_PTR(m_omega22_poly[ipoly]));
m_astar(i,j) = poly8(z, DATA_PTR(m_astar_poly[ipoly]));
m_bstar(i,j) = poly8(z, DATA_PTR(m_bstar_poly[ipoly]));
m_cstar(i,j) = poly8(z, DATA_PTR(m_cstar_poly[ipoly]));

14
src/transport/TransportFactory.cpp
View file

@ -214,34 +214,34 @@ Transport* TransportFactory::newTransport(const std::string& transportModel,
break;
case cMulticomponent:
tr = new MultiTransport;
dynamic_cast<GasTransport*>(tr)->init(phase, 0, log_level);
tr->init(phase, 0, log_level);
break;
case CK_Multicomponent:
tr = new MultiTransport;
dynamic_cast<GasTransport*>(tr)->init(phase, CK_Mode, log_level);
tr->init(phase, CK_Mode, log_level);
break;
case cMixtureAveraged:
tr = new MixTransport;
dynamic_cast<GasTransport*>(tr)->init(phase, 0, log_level);
tr->init(phase, 0, log_level);
break;
case CK_MixtureAveraged:
tr = new MixTransport;
dynamic_cast<GasTransport*>(tr)->init(phase, CK_Mode, log_level);
tr->init(phase, CK_Mode, log_level);
break;
case cHighP:
tr = new HighPressureGasTransport;
dynamic_cast<GasTransport*>(tr)->init(phase, 0, log_level);
tr->init(phase, 0, log_level);
break;
case cSolidTransport:
tr = new SolidTransport;
initSolidTransport(tr, phase, log_level);
dynamic_cast<GasTransport*>(tr)->setThermo(*phase);
tr->setThermo(*phase);
break;
case cDustyGasTransport:
tr = new DustyGasTransport;
gastr = new MultiTransport;
dynamic_cast<GasTransport*>(gastr)->init(phase, 0, log_level);
gastr->init(phase, 0, log_level);
dtr = (DustyGasTransport*)tr;
dtr->initialize(phase, gastr);
break;

2
src/transport/TransportParams.cpp
View file

@ -51,6 +51,8 @@ GasTransportParams::GasTransportParams() :
dipole(0, 0),
delta(0, 0)
{
warn_deprecated("class GasTransportParams",
"To be removed after Cantera 2.2.");
}
} // End of namespace Cantera