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

Merged viscosity calculations from MixTransport and MultiTransport

Browse source
This commit is contained in:
Ray Speth 2012-05-11 15:27:44 +00:00
parent 74e9fa050b
commit f8308d5853
6 changed files with 394 additions and 618 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

161
include/cantera/transport/GasTransport.h Normal file
View file

@ -0,0 +1,161 @@
/**
* @file GasTransport.h
*/
#ifndef CT_GAS_TRANSPORT_H
#define CT_GAS_TRANSPORT_H
#include "TransportBase.h"
namespace Cantera {
//! Class GasTransport implements some functions and properties that are
//! shared by the MixTransport and MultiTransport classes.
class GasTransport : public Transport
{
public:
virtual ~GasTransport() {}
GasTransport(const GasTransport& right);
GasTransport& operator=(const GasTransport& right);
//! Viscosity of the mixture (kg /m /s)
/*!
* The viscosity is computed using the Wilke mixture rule (kg /m /s)
*
* \f[
* \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}.
* \f]
*
* Here \f$ \mu_k \f$ is the viscosity of pure species \e k, and
*
* \f[
* \Phi_{k,j} = \frac{\left[1
* + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2}
* {\sqrt{8}\sqrt{1 + M_k/M_j}}
* \f]
*
* @return Returns the viscosity of the mixture ( units = Pa s = kg /m /s)
*
* @see updateViscosity_T();
*/
virtual doublereal viscosity();
//! Get the pure-species viscosities
virtual void getSpeciesViscosities(doublereal* const visc) {
update_T();
updateViscosity_T();
std::copy(m_visc.begin(), m_visc.end(), visc);
}
protected:
GasTransport(ThermoPhase* thermo=0);
virtual bool initGas(GasTransportParams& tr);
virtual void update_T();
virtual void update_C() = 0;
//! Update the temperature-dependent viscosity terms.
/**
* Updates the array of pure species viscosities, and the weighting
* functions in the viscosity mixture rule. The flag m_visc_ok is set to true.
*
* The formula for the weighting function is from Poling and Prausnitz,
* Eq. (9-5.14):
* \f[
* \phi_{ij} = \frac{ \left[ 1 + \left( \mu_i / \mu_j \right)^{1/2} \left( M_j / M_i \right)^{1/4} \right]^2 }
* {\left[ 8 \left( 1 + M_i / M_j \right) \right]^{1/2}}
* \f]
*/
virtual void updateViscosity_T();
//! Update the pure-species viscosities. These are evaluated from the
//! polynomial fits of the temperature and are assumed to be independent
//! of pressure.
virtual void updateSpeciesViscosities();
//! Vector of species mole fractions. These are processed so that all mole
//! fractions are >= MIN_X. Length = m_kk.
vector_fp m_molefracs;
//! Internal storage for the viscosity of the mixture (kg /m /s)
doublereal m_viscmix;
//! Update boolean for mixture rule for the mixture viscosity
bool m_visc_ok;
//! Update boolean for the weighting factors for the mixture viscosity
bool m_viscwt_ok;
//! Update boolean for the species viscosities
bool m_spvisc_ok;
//! Type of the polynomial fits to temperature. CK_Mode means Chemkin mode.
//! Currently CA_Mode is used which are different types of fits to temperature.
int m_mode;
//! m_phi is a Viscosity Weighting Function. size = m_nsp * n_nsp
DenseMatrix m_phi;
//! work space length = m_kk
vector_fp m_spwork;
//! vector of species viscosities (kg /m /s). These are used in Wilke's
//! rule to calculate the viscosity of the solution. length = m_kk.
vector_fp m_visc;
//! Polynomial fits to the viscosity of each species. m_visccoeffs[k] is
//! the vector of polynomial coefficients for species k that fits the
//! viscosity as a function of temperature.
std::vector<vector_fp> m_visccoeffs;
//! Local copy of the species molecular weights.
vector_fp m_mw;
//! Holds square roots of molecular weight ratios
/*!
* m_wratjk(j,k) = sqrt(mw[j]/mw[k]) j < k
* m_wratjk(k,j) = sqrt(sqrt(mw[j]/mw[k])) j < k
*/
DenseMatrix m_wratjk;
//! Holds square roots of molecular weight ratios
/*!
* m_wratjk1(j,k) = sqrt(1.0 + mw[k]/mw[j]) j < k
*/
DenseMatrix m_wratkj1;
//! vector of square root of species viscosities sqrt(kg /m /s). These are
//! used in Wilke's rule to calculate the viscosity of the solution.
//! length = m_kk.
vector_fp m_sqvisc;
//! Powers of the ln temperature, up to fourth order
vector_fp m_polytempvec;
//! Current value of the temperature at which the properties in this object
//! are calculated (Kelvin).
doublereal m_temp;
//! Current value of Boltzman's constant times the temperature (Joules)
doublereal m_kbt;
//! current value of Boltzman's constant times the temperature.
//! (Joules) to 1/2 power
doublereal m_sqrt_kbt;
//! current value of temperature to 1/2 power
doublereal m_sqrt_t;
//! Current value of the log of the temperature
doublereal m_logt;
//! Current value of temperature to 1/4 power
doublereal m_t14;
//! Current value of temperature to the 3/2 power
doublereal m_t32;
};
} // namespace Cantera
#endif

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

@ -18,16 +18,14 @@
#include <algorithm>
// Cantera includes
#include "TransportBase.h"
#include "GasTransport.h"
#include "cantera/numerics/DenseMatrix.h"
namespace Cantera
{
class GasTransportParams;
//! Class MixTransport implements mixture-averaged transport properties for ideal gas mixtures.
/*!
* The model is based on that described by Kee, Coltrin, and Glarborg, "Theoretical and
@ -66,15 +64,12 @@ class GasTransportParams;
*
*
*/
class MixTransport : public Transport
class MixTransport : public GasTransport
{
protected:
//! Default constructor.
/*!
*
*/
MixTransport();
public:
@ -106,9 +101,8 @@ public:
*/
virtual Transport* duplMyselfAsTransport() const;
//! Destructor
virtual ~MixTransport();
virtual ~MixTransport() {}
//! Return the model id for transport
/*!
@ -118,38 +112,6 @@ public:
return cMixtureAveraged;
}
//! Viscosity of the mixture (kg /m /s)
/*!
* The viscosity is computed using the Wilke mixture rule (kg /m /s)
*
* \f[
* \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}.
* \f]
*
* Here \f$ \mu_k \f$ is the viscosity of pure species \e k, and
*
* \f[
* \Phi_{k,j} = \frac{\left[1
* + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2}
* {\sqrt{8}\sqrt{1 + M_k/M_j}}
* \f]
*
* @return Returns the viscosity of the mixture ( units = Pa s = kg /m /s)
*
* @see updateViscosity_T();
*/
virtual doublereal viscosity();
//! returns the vector of species viscosities
/*!
* @param visc Vector of species viscosities
*/
virtual void getSpeciesViscosities(doublereal* const visc) {
update_T();
updateViscosity_T();
copy(m_visc.begin(), m_visc.end(), visc);
}
//! Return the thermal diffusion coefficients
/*!
* For this approximation, these are all zero.
@ -285,7 +247,6 @@ public:
friend class TransportFactory;
//! Return a structure containing all of the pertinent parameters about a species that was
//! used to construct the Transport properties in this object.
/*!
@ -295,8 +256,6 @@ public:
*/
struct GasTransportData getGasTransportData(int kspec) const;
private:
//! Calculate the pressure from the ideal gas law
@ -305,22 +264,6 @@ private:
m_thermo->temperature());
}
//! Update the temperature-dependent viscosity terms.
/*!
* Updates the array of pure species viscosities, and the weighting functions in the viscosity mixture rule.
* The flag m_visc_ok is set to true.
*
* The formula for the weighting function is from Poling and Prausnitz.
* See Eq. (9-5.14) of Poling, Prausnitz, and O'Connell. The equation for the weighting function
* \f$ \phi_{ij} \f$ is reproduced below.
*
* \f[
* \phi_{ij} = \frac{ \left[ 1 + \left( \mu_i / \mu_j \right)^{1/2} \left( M_j / M_i \right)^{1/4} \right]^2 }
* {\left[ 8 \left( 1 + M_i / M_j \right) \right]^{1/2}}
* \f]
*/
void updateViscosity_T();
//! Update the temperature dependent parts of the species thermal conductivities
/*!
* These are evaluated from the polynomial fits of the temperature and are assumed to be
@ -328,13 +271,6 @@ private:
*/
void updateCond_T();
//! Update the species viscosities
/*!
* These are evaluated from the polynomial fits of the temperature and are assumed to be
* independent of pressure
*/
void updateSpeciesViscosities();
//! Update the binary diffusion coefficients
/*!
* These are evaluated from the polynomial fits of the temperature at the unit pressure of 1 Pa.
@ -344,21 +280,6 @@ private:
// --------- Member Data -------------
private:
//! Minimum value of the temperature that this transport parameterization is valid
doublereal m_tmin;
//! Maximum value of the temperature that this transport parameterization is valid
doublereal m_tmax;
//! Local copy of the species molecular weights.
vector_fp m_mw;
//! Polynomial fits to the viscosity of each species
/*!
* m_visccoeffs[k] is vector of polynomial coefficients for species k
* that fits the viscosity as a function of temperature
*/
std::vector<vector_fp> m_visccoeffs;
//! Polynomial fits to the thermal conductivity of each species
/*!
@ -383,32 +304,12 @@ private:
*/
std::vector<vector_fp> m_diffcoeffs;
//! Powers of the ln temperature
/*!
* up to fourth order
*/
vector_fp m_polytempvec;
//! Matrix of binary diffusion coefficients at the reference pressure and the current temperature
/*!
* Size is nsp x nsp
*/
DenseMatrix m_bdiff;
//! vector of species viscosities (kg /m /s)
/*!
* These are used in wilke's rule to calculate the viscosity of the solution
* length = m_kk
*/
vector_fp m_visc;
//! vector of square root of species viscosities sqrt(kg /m /s)
/*!
* These are used in wilke's rule to calculate the viscosity of the solution
* length = m_kk
*/
vector_fp m_sqvisc;
//! vector of species thermal conductivities (W/m /K)
/*!
* These are used in wilke's rule to calculate the viscosity of the solution
@ -417,78 +318,12 @@ private:
*/
vector_fp m_cond;
//! Vector of species molefractions
/*!
* These are processed so that all mole fractions are >= MIN_X
* Length = m_kk
*/
vector_fp m_molefracs;
//! m_phi is a Viscosity Weighting Function
/*!
* size = m_nsp * n_nsp
*/
DenseMatrix m_phi;
//! Holds square roots or molecular weight ratios
/*!
* m_wratjk(j,k) = sqrt(mw[j]/mw[k]) j < k
* m_wratjk(k,j) = sqrt(sqrt(mw[j]/mw[k])) j < k
*/
DenseMatrix m_wratjk;
//! Holds square roots of molecular weight ratios
/*!
* m_wratjk1(j,k) = sqrt(1.0 + mw[k]/mw[j]) j < k
*/
DenseMatrix m_wratkj1;
//! Current value of the temperature at which the properties in this object are calculated (Kelvin)
doublereal m_temp;
//! Current value of the log of the temperature
doublereal m_logt;
//! Current value of Boltzman's constant times the temperature (Joules)
doublereal m_kbt;
//! Current value of temperature to 1/4 power
doublereal m_t14;
//! Current value of temperature to the 3/2 power
doublereal m_t32;
//! current value of Boltzman's constant times the temperature (Joules) to 1/2 power
doublereal m_sqrt_kbt;
//! current value of temperature to 1/2 power
doublereal m_sqrt_t;
//! Internal storage for the calculated mixture thermal conductivity
/*!
* Units = W /m /K
*/
doublereal m_lambda;
//! Internal storage for the viscosity of the mixture (kg /m /s)
doublereal m_viscmix;
//! work space length = m_kk
vector_fp m_spwork;
//! Update boolean for mixture rule for the mixture viscosity
bool m_viscmix_ok;
//! Update boolean for the weighting factors for the mixture viscosity
bool m_viscwt_ok;
//! Update boolean for the species viscosities
bool m_spvisc_ok;
//! Update boolean for the binary diffusivities at unit pressure
bool m_bindiff_ok;
@ -498,13 +333,6 @@ private:
//! Update boolean for the mixture rule for the mixture thermal conductivity
bool m_condmix_ok;
//! Type of the polynomial fits to temperature
/*!
* CK_Mode means chemkin mode. Currently CA_Mode is used which are different types
* of fits to temperature.
*/
int m_mode;
//! Lennard-Jones well-depth of the species in the current phase
/*!
* Not used in this routine -> just a passthrough

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

@ -13,7 +13,7 @@
//#undef CHEMKIN_COMPATIBILITY_MODE
// Cantera includes
#include "TransportBase.h"
#include "GasTransport.h"
#include "cantera/numerics/DenseMatrix.h"
namespace Cantera
@ -78,7 +78,7 @@ public:
*
* @ingroup transportProps
*/
class MultiTransport : public Transport
class MultiTransport : public GasTransport
{
protected:
@ -92,7 +92,7 @@ protected:
public:
//! Destructor
virtual ~MultiTransport();
virtual ~MultiTransport() {}
// overloaded base class methods
virtual int model() const {
@ -103,14 +103,6 @@ public:
}
}
virtual doublereal viscosity();
virtual void getSpeciesViscosities(doublereal* const visc) {
updateViscosity_T();
std::copy(m_visc.begin(), m_visc.end(), visc);
}
//! Return the thermal diffusion coefficients (kg/m/s)
/*!
* Eqn. (12.126) displays how they are calculated. The reference work is from
@ -211,10 +203,6 @@ public:
*/
DEPRECATED(virtual void setOptions_GMRES(int m, doublereal eps));
/**
* @internal
*/
//! Initialize the transport operator with parameters from GasTransportParams object
/*!
* @param tr input GasTransportParams object
@ -233,24 +221,15 @@ public:
protected:
//! Update basic temperature-dependent quantities if the temperature has changed.
void updateTransport_T();
void update_T();
//! Update basic concentration-dependent quantities if the concentrations have changed.
void updateTransport_C();
void update_C();
//! Update the temperature-dependent terms needed to compute the thermal
//! conductivity and thermal diffusion coefficients.
void updateThermal_T();
//! Update the temperature-dependent viscosity terms
void updateViscosity_T();
//! Update the temperature-dependent viscosity terms.
//! Updates the array of pure species viscosities and the weighting
//! functions in the viscosity mixture rule.
//! The flag m_visc_ok is set to true.
void updateSpeciesViscosities_T();
//! Update the binary diffusion coefficients.
//! These are evaluated from the polynomial fits at unit pressure (1 Pa).
void updateDiff_T();
@ -258,26 +237,16 @@ protected:
private:
doublereal m_diff_tlast;
doublereal m_spvisc_tlast;
doublereal m_visc_tlast;
doublereal m_thermal_tlast;
doublereal m_tmin;
doublereal m_tmax;
vector_fp m_mw;
// polynomial fits
std::vector<vector_fp> m_visccoeffs;
std::vector<vector_fp> m_diffcoeffs;
vector_fp m_polytempvec;
std::vector<vector_fp> m_diffcoeffs;
// property values
DenseMatrix m_bdiff;
vector_fp m_visc;
vector_fp m_sqvisc;
vector_fp m_molefracs;
std::vector<std::vector<int> > m_poly;
std::vector<vector_fp> m_astar_poly;
@ -297,9 +266,6 @@ private:
//! Dense matrix for omega22
DenseMatrix m_om22;
DenseMatrix m_phi; // viscosity weighting functions
DenseMatrix m_wratjk, m_wratkj1;
vector_fp m_zrot;
vector_fp m_crot;
vector_fp m_cinternal;
@ -307,9 +273,6 @@ private:
vector_fp m_alpha;
vector_fp m_dipoleDiag;
doublereal m_temp, m_logt, m_kbt, m_t14, m_t32;
doublereal m_sqrt_kbt, m_sqrt_t;
vector_fp m_sqrt_eps_k;
DenseMatrix m_log_eps_k;
vector_fp m_frot_298;
@ -329,18 +292,15 @@ private:
doublereal m_eps_gmres; //!< @deprecated
// work space
vector_fp m_spwork, m_spwork1, m_spwork2, m_spwork3;
vector_fp m_spwork1, m_spwork2, m_spwork3;
void correctBinDiffCoeffs();
//! Boolean indicating viscosity is up to date
bool m_visc_ok;
bool m_spvisc_ok;
bool m_diff_ok;
bool m_abc_ok;
bool m_l0000_ok;
bool m_lmatrix_soln_ok;
int m_mode;
//! Evaluate the L0000 matrices
/*!
@ -356,9 +316,6 @@ private:
void eval_L0010(const doublereal* const x);
//! Evaluate the L1000 matrices
/*!
*
*/
void eval_L1000();
void eval_L0100();

206
src/transport/GasTransport.cpp Normal file
View file

@ -0,0 +1,206 @@
#include "cantera/transport/GasTransport.h"
#include "cantera/transport/TransportParams.h"
namespace Cantera {
GasTransport::GasTransport(ThermoPhase* thermo) :
Transport(thermo),
m_molefracs(0),
m_viscmix(0.0),
m_visc_ok(false),
m_viscwt_ok(false),
m_spvisc_ok(false),
m_mode(0),
m_phi(0,0),
m_spwork(0),
m_visc(0),
m_visccoeffs(0),
m_mw(0),
m_wratjk(0,0),
m_wratkj1(0,0),
m_sqvisc(0),
m_polytempvec(5),
m_temp(-1.0),
m_kbt(0.0),
m_sqrt_kbt(0.0),
m_sqrt_t(0.0),
m_logt(0.0),
m_t14(0.0),
m_t32(0.0)
{
}
GasTransport::GasTransport(const GasTransport& right) :
m_molefracs(0),
m_viscmix(0.0),
m_visc_ok(false),
m_viscwt_ok(false),
m_spvisc_ok(false),
m_mode(0),
m_phi(0,0),
m_spwork(0),
m_visc(0),
m_visccoeffs(0),
m_mw(0),
m_wratjk(0,0),
m_wratkj1(0,0),
m_sqvisc(0),
m_polytempvec(5),
m_temp(-1.0),
m_kbt(0.0),
m_sqrt_kbt(0.0),
m_sqrt_t(0.0),
m_logt(0.0),
m_t14(0.0),
m_t32(0.0)
{
}
GasTransport& GasTransport::operator=(const GasTransport& right)
{
m_molefracs = right.m_molefracs;
m_viscmix = right.m_viscmix;
m_visc_ok = right.m_visc_ok;
m_viscwt_ok = right.m_viscwt_ok;
m_spvisc_ok = right.m_spvisc_ok;
m_mode = right.m_mode;
m_phi = right.m_phi;
m_spwork = right.m_spwork;
m_visc = right.m_visc;
m_mw = right.m_mw;
m_wratjk = right.m_wratjk;
m_wratkj1 = right.m_wratkj1;
m_sqvisc = right.m_sqvisc;
m_polytempvec = right.m_polytempvec;
m_temp = right.m_temp;
m_kbt = right.m_kbt;
m_sqrt_kbt = right.m_sqrt_kbt;
m_sqrt_t = right.m_sqrt_t;
m_logt = right.m_logt;
m_t14 = right.m_t14;
m_t32 = right.m_t32;
return *this;
}
bool GasTransport::initGas(GasTransportParams& tr)
{
// constant mixture attributes
m_thermo = tr.thermo;
m_nsp = m_thermo->nSpecies();
m_molefracs.resize(m_nsp);
m_spwork.resize(m_nsp);
m_visc.resize(m_nsp);
m_phi.resize(m_nsp, m_nsp, 0.0);
// 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;
}
void GasTransport::update_T(void) {
m_temp = m_thermo->temperature();
m_kbt = Boltzmann * m_temp;
m_sqrt_kbt = sqrt(Boltzmann*m_temp);
m_logt = log(m_temp);
m_sqrt_t = sqrt(m_temp);
m_t14 = sqrt(m_sqrt_t);
m_t32 = m_temp * m_sqrt_t;
// compute powers of log(T)
m_polytempvec[0] = 1.0;
m_polytempvec[1] = m_logt;
m_polytempvec[2] = m_logt*m_logt;
m_polytempvec[3] = m_logt*m_logt*m_logt;
m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt;
// temperature has changed, so polynomial fits will need to be redone
m_visc_ok = false;
m_spvisc_ok = false;
m_viscwt_ok = false;
}
doublereal GasTransport::viscosity()
{
update_T();
update_C();
if (m_visc_ok) {
return m_viscmix;
}
doublereal vismix = 0.0;
// update m_visc and m_phi if necessary
if (!m_viscwt_ok) {
updateViscosity_T();
}
multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork));
for (size_t k = 0; k < m_nsp; k++) {
vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
}
m_viscmix = vismix;
return vismix;
}
void GasTransport::updateViscosity_T()
{
doublereal vratiokj, wratiojk, factor1;
if (!m_spvisc_ok) {
updateSpeciesViscosities();
}
// see Eq. (9-5.15) of Reid, Prausnitz, and Poling
for (size_t j = 0; j < m_nsp; j++) {
for (size_t k = j; k < m_nsp; k++) {
vratiokj = m_visc[k]/m_visc[j];
wratiojk = m_mw[j]/m_mw[k];
// Note that m_wratjk(k,j) holds the square root of m_wratjk(j,k)!
factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
m_phi(k,j) = factor1*factor1 / (SqrtEight * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
}
m_viscwt_ok = true;
}
void GasTransport::updateSpeciesViscosities()
{
if (m_mode == CK_Mode) {
for (size_t k = 0; k < m_nsp; k++) {
m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
m_sqvisc[k] = sqrt(m_visc[k]);
}
} else {
for (size_t k = 0; k < m_nsp; k++) {
// the polynomial fit is done for sqrt(visc/sqrt(T))
m_sqvisc[k] = m_t14 * dot5(m_polytempvec, m_visccoeffs[k]);
m_visc[k] = (m_sqvisc[k] * m_sqvisc[k]);
}
}
m_spvisc_ok = true;
}
}

230
src/transport/MixTransport.cpp
View file

@ -29,38 +29,14 @@ namespace Cantera
//====================================================================================================================
MixTransport::MixTransport() :
m_tmin(-1.0),
m_tmax(100000.),
m_mw(0),
m_visccoeffs(0),
m_condcoeffs(0),
m_diffcoeffs(0),
m_polytempvec(0),
m_bdiff(0, 0),
m_visc(0),
m_sqvisc(0),
m_cond(0),
m_molefracs(0),
m_phi(0,0),
m_wratjk(0,0),
m_wratkj1(0,0),
m_temp(-1.0),
m_logt(0.0),
m_kbt(0.0),
m_t14(0.0),
m_t32(0.0),
m_sqrt_kbt(0.0),
m_sqrt_t(0.0),
m_lambda(0.0),
m_viscmix(0.0),
m_spwork(0),
m_viscmix_ok(false),
m_viscwt_ok(false),
m_spvisc_ok(false),
m_bindiff_ok(false),
m_spcond_ok(false),
m_condmix_ok(false),
m_mode(0),
m_eps(0),
m_diam(0, 0),
m_dipoleDiag(0),
@ -72,38 +48,15 @@ MixTransport::MixTransport() :
}
//====================================================================================================================
MixTransport::MixTransport(const MixTransport& right) :
m_tmin(-1.0),
m_tmax(100000.),
m_mw(0),
m_visccoeffs(0),
GasTransport(right),
m_condcoeffs(0),
m_diffcoeffs(0),
m_polytempvec(0),
m_bdiff(0, 0),
m_visc(0),
m_sqvisc(0),
m_cond(0),
m_molefracs(0),
m_phi(0,0),
m_wratjk(0,0),
m_wratkj1(0,0),
m_temp(-1.0),
m_logt(0.0),
m_kbt(0.0),
m_t14(0.0),
m_t32(0.0),
m_sqrt_kbt(0.0),
m_sqrt_t(0.0),
m_lambda(0.0),
m_viscmix(0.0),
m_spwork(0),
m_viscmix_ok(false),
m_viscwt_ok(false),
m_spvisc_ok(false),
m_bindiff_ok(false),
m_spcond_ok(false),
m_condmix_ok(false),
m_mode(0),
m_eps(0),
m_diam(0, 0),
m_dipoleDiag(0),
@ -127,40 +80,16 @@ MixTransport& MixTransport::operator=(const MixTransport& right)
if (&right == this) {
return *this;
}
Transport::operator=(right);
GasTransport::operator=(right);
m_tmin = right.m_tmin;
m_tmax = right.m_tmax;
m_mw =right.m_mw;
m_visccoeffs = right.m_visccoeffs;
m_condcoeffs = right.m_condcoeffs;
m_diffcoeffs = right.m_diffcoeffs;
m_polytempvec = right.m_polytempvec;
m_bdiff = right.m_bdiff;
m_visc = right.m_visc;
m_sqvisc = right.m_sqvisc;
m_cond = right.m_cond;
m_molefracs = right.m_molefracs;
m_phi = right.m_phi;
m_wratjk = right.m_wratjk;
m_wratkj1 = right.m_wratkj1;
m_temp = right.m_temp;
m_logt = right.m_logt;
m_kbt = right.m_kbt;
m_t14 = right.m_t14;
m_t32 = right.m_t32;
m_sqrt_kbt = right.m_sqrt_kbt;
m_sqrt_t = right.m_sqrt_t;
m_lambda = right.m_lambda;
m_viscmix = right.m_viscmix;
m_spwork = right.m_spwork;
m_viscmix_ok = right.m_viscmix_ok;
m_viscwt_ok = right.m_viscwt_ok;
m_spvisc_ok = right.m_spvisc_ok;
m_bindiff_ok = right.m_bindiff_ok;
m_spcond_ok = right.m_spcond_ok;
m_condmix_ok = right.m_condmix_ok;
m_mode = right.m_mode;
m_eps = right.m_eps;
m_diam = right.m_diam;
m_dipoleDiag = right.m_dipoleDiag;
@ -186,24 +115,11 @@ Transport* MixTransport::duplMyselfAsTransport() const
MixTransport* tr = new MixTransport(*this);
return (dynamic_cast<Transport*>(tr));
}
//====================================================================================================================
MixTransport::~MixTransport()
{
}
//====================================================================================================================
bool MixTransport::initGas(GasTransportParams& tr)
{
// constant substance attributes
m_thermo = tr.thermo;
m_nsp = m_thermo->nSpecies();
m_tmin = m_thermo->minTemp();
m_tmax = m_thermo->maxTemp();
// make a local copy of the molecular weights
m_mw.resize(m_nsp);
copy(m_thermo->molecularWeights().begin(),
m_thermo->molecularWeights().end(), m_mw.begin());
GasTransport::initGas(tr);
// copy polynomials and parameters into local storage
m_visccoeffs = tr.visccoeffs;
@ -221,78 +137,16 @@ bool MixTransport::initGas(GasTransportParams& tr)
m_dipoleDiag[i] = tr.dipole(i,i);
}
m_phi.resize(m_nsp, m_nsp, 0.0);
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
size_t j, k;
for (j = 0; j < m_nsp; j++)
for (k = j; k < m_nsp; k++) {
m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
m_wratjk(k,j) = sqrt(m_wratjk(j,k));
m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
}
m_polytempvec.resize(5);
m_visc.resize(m_nsp);
m_sqvisc.resize(m_nsp);
m_cond.resize(m_nsp);
m_bdiff.resize(m_nsp, m_nsp);
m_molefracs.resize(m_nsp);
m_spwork.resize(m_nsp);
// set flags all false
m_viscmix_ok = false;
m_viscwt_ok = false;
m_spvisc_ok = false;
m_spcond_ok = false;
m_condmix_ok = false;
return true;
}
//====================================================================================================================
// Viscosity of the mixture
/*
*
* The viscosity is computed using the Wilke mixture rule.
*
* \f[
* \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}.
* \f]
*
* Here \f$ \mu_k \f$ is the viscosity of pure species \e k, and
*
* \f[
* \Phi_{k,j} = \frac{\left[1
* + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2}
* {\sqrt{8}\sqrt{1 + M_k/M_j}}
* \f]
*
* @see updateViscosity_T();
*/
doublereal MixTransport::viscosity()
{
update_T();
update_C();
if (m_viscmix_ok) {
return m_viscmix;
}
doublereal vismix = 0.0;
// update m_visc and m_phi if necessary
if (!m_viscwt_ok) {
updateViscosity_T();
}
multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork));
for (size_t k = 0; k < m_nsp; k++) {
vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
}
m_viscmix = vismix;
return vismix;
}
//====================================================================================================================
// Returns the matrix of binary diffusion coefficients.
/*
@ -498,26 +352,8 @@ void MixTransport::update_T()
throw CanteraError("MixTransport::update_T",
"negative temperature "+fp2str(t));
}
m_temp = t;
m_logt = log(m_temp);
m_kbt = Boltzmann * m_temp;
m_sqrt_t = sqrt(m_temp);
m_t14 = sqrt(m_sqrt_t);
m_t32 = m_temp * m_sqrt_t;
m_sqrt_kbt = sqrt(Boltzmann*m_temp);
// compute powers of log(T)
m_polytempvec[0] = 1.0;
m_polytempvec[1] = m_logt;
m_polytempvec[2] = m_logt*m_logt;
m_polytempvec[3] = m_logt*m_logt*m_logt;
m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt;
// temperature has changed, so polynomial fits will need to be
// redone.
m_viscmix_ok = false;
m_spvisc_ok = false;
m_viscwt_ok = false;
GasTransport::update_T();
// temperature has changed, so polynomial fits will need to be redone.
m_spcond_ok = false;
m_bindiff_ok = false;
m_condmix_ok = false;
@ -534,7 +370,7 @@ void MixTransport::update_C()
// be recomputed before use, and update the local mole
// fractions.
m_viscmix_ok = false;
m_visc_ok = false;
m_condmix_ok = false;
m_thermo->getMoleFractions(DATA_PTR(m_molefracs));
@ -597,60 +433,8 @@ void MixTransport::updateDiff_T()
/*
* Update the pure-species viscosities.
*/
void MixTransport::updateSpeciesViscosities()
{
if (m_mode == CK_Mode) {
for (size_t k = 0; k < m_nsp; k++) {
m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
m_sqvisc[k] = sqrt(m_visc[k]);
}
} else {
for (size_t k = 0; k < m_nsp; k++) {
// the polynomial fit is done for sqrt(visc/sqrt(T))
m_sqvisc[k] = m_t14 * dot5(m_polytempvec, m_visccoeffs[k]);
m_visc[k] = (m_sqvisc[k] * m_sqvisc[k]);
}
}
m_spvisc_ok = true;
}
//====================================================================================================================
// Update the temperature-dependent viscosity terms.
/*
* Updates the array of pure species viscosities, and the weighting functions in the viscosity mixture rule.
* The flag m_visc_ok is set to true.
*
* The formula for the weighting function is from Poling and Prausnitz.
* See Eq. (9-5.14) of Poling, Prausnitz, and O'Connell. The equation for the weighting function
* \f$ \phi_{ij} \f$ is reproduced below.
*
* \f[
* \phi_{ij} = \frac{ \left[ 1 + \left( \mu_i / \mu_j \right)^{1/2} \left( M_j / M_i \right)^{1/4} \right]^2 }
* {\left[ 8 \left( 1 + M_i / M_j \right) \right]^{1/2}}
* \f]
*/
void MixTransport::updateViscosity_T()
{
doublereal vratiokj, wratiojk, factor1;
if (!m_spvisc_ok) {
updateSpeciesViscosities();
}
// see Eq. (9-5.15) of Reid, Prausnitz, and Poling
for (size_t j = 0; j < m_nsp; j++) {
for (size_t k = j; k < m_nsp; k++) {
vratiokj = m_visc[k]/m_visc[j];
wratiojk = m_mw[j]/m_mw[k];
// Note that m_wratjk(k,j) holds the square root of
// m_wratjk(j,k)!
factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
m_phi(k,j) = factor1*factor1 / (SqrtEight * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
}
m_viscwt_ok = true;
}
//====================================================================================================================
/*
* This function returns a Transport data object for a given species.

180
src/transport/MultiTransport.cpp
View file

@ -29,23 +29,13 @@ using namespace std;
* Mole fractions below MIN_X will be set to MIN_X when computing
* transport properties.
*/
#define MIN_X 1.e-20
namespace Cantera
{
/////////////////////////// constants //////////////////////////
// const doublereal ThreeSixteenths = 3.0/16.0;
///////////////////// helper functions /////////////////////////
/**
* @internal
*
@ -92,33 +82,15 @@ void L_Matrix::mult(const doublereal* b, doublereal* prod) const
//////////////////// class MultiTransport methods //////////////
MultiTransport::MultiTransport(thermo_t* thermo)
: Transport(thermo),
m_temp(-1.0)
: GasTransport(thermo)
{
}
MultiTransport::~MultiTransport()
{
}
//====================================================================================================================
bool MultiTransport::initGas(GasTransportParams& tr)
{
// constant mixture attributes
//m_phase = tr.mix;
m_thermo = tr.thermo;
m_nsp = m_thermo->nSpecies();
m_tmin = m_thermo->minTemp();
m_tmax = m_thermo->maxTemp();
// make a local copy of the molecular weights
m_mw.resize(m_nsp);
copy(m_thermo->molecularWeights().begin(),
m_thermo->molecularWeights().end(), m_mw.begin());
GasTransport::initGas(tr);
// copy polynomials and parameters into local storage
m_poly = tr.poly;
@ -149,21 +121,8 @@ bool MultiTransport::initGas(GasTransportParams& tr)
m_frot_298.resize(m_nsp);
m_rotrelax.resize(m_nsp);
m_phi.resize(m_nsp, m_nsp, 0.0);
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
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_cinternal.resize(m_nsp);
m_polytempvec.resize(5);
m_visc.resize(m_nsp);
m_sqvisc.resize(m_nsp);
m_bdiff.resize(m_nsp, m_nsp);
//m_poly.resize(m_nsp);
@ -172,19 +131,13 @@ bool MultiTransport::initGas(GasTransportParams& tr)
m_bstar.resize(m_nsp, m_nsp);
m_cstar.resize(m_nsp, m_nsp);
m_molefracs.resize(m_nsp);
// set flags all false
m_visc_ok = false;
m_spvisc_ok = false;
m_diff_ok = false;
m_abc_ok = false;
m_l0000_ok = false;
m_lmatrix_soln_ok = false;
m_diff_tlast = 0.0;
m_spvisc_tlast = 0.0;
m_visc_tlast = 0.0;
m_thermal_tlast = 0.0;
// use LU decomposition by default
@ -195,12 +148,10 @@ bool MultiTransport::initGas(GasTransportParams& tr)
m_eps_gmres = 1.e-4;
// some work space
m_spwork.resize(m_nsp);
m_spwork1.resize(m_nsp);
m_spwork2.resize(m_nsp);
m_spwork3.resize(m_nsp);
// precompute and store log(epsilon_ij/k_B)
m_log_eps_k.resize(m_nsp, m_nsp);
// int j;
@ -211,7 +162,6 @@ bool MultiTransport::initGas(GasTransportParams& tr)
}
}
// precompute and store constant parts of the Parker rotational
// collision number temperature correction
const doublereal sq298 = sqrt(298.0);
@ -223,46 +173,9 @@ bool MultiTransport::initGas(GasTransportParams& tr)
m_sqrt_eps_k[k]/sq298);
}
// // install updaters
// m_update_transport_T = m_thermo->installUpdater_T(
// new UpdateTransport_T<MultiTransport>(*this));
// m_update_transport_C = m_thermo->installUpdater_C(
// new UpdateTransport_C<MultiTransport>(*this));
// m_update_spvisc_T = m_thermo->installUpdater_T(
// new UpdateSpeciesVisc<MultiTransport>(*this));
// m_update_visc_T = m_thermo->installUpdater_T(
// new UpdateVisc_T<MultiTransport>(*this));
// m_update_diff_T = m_thermo->installUpdater_T(
// new UpdateDiff_T<MultiTransport>(*this));
// m_update_thermal_T = m_thermo->installUpdater_T(
// new UpdateThermal_T<MultiTransport>(*this));
return true;
}
/****************** viscosity ******************************/
doublereal MultiTransport::viscosity()
{
doublereal vismix = 0.0, denom;
// update m_visc if necessary
updateViscosity_T();
// update the mole fractions
updateTransport_C();
for (size_t k = 0; k < m_nsp; k++) {
denom = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
denom += m_phi(k,j) * m_molefracs[j];
}
vismix += m_molefracs[k] * m_visc[k]/denom;
}
return vismix;
}
//====================================================================================================================
/******************* binary diffusion coefficients **************/
@ -282,7 +195,6 @@ void MultiTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d)
}
/****************** thermal conductivity **********************/
/**
@ -324,7 +236,7 @@ void MultiTransport::solveLMatrixEquation()
// properties.
updateThermal_T();
updateTransport_C();
update_C();
// Copy the mole fractions twice into the last two blocks of
// the right-hand-side vector m_b. The first block of m_b was
@ -700,7 +612,7 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
doublereal p = pressure_ig();
// update the mole fractions
updateTransport_C();
update_C();
// update the binary diffusion coefficients
updateDiff_T();
@ -736,9 +648,8 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
void MultiTransport::getMixDiffCoeffs(doublereal* const d)
{
// update the mole fractions
updateTransport_C();
update_C();
// update the binary diffusion coefficients if necessary
updateDiff_T();
@ -768,39 +679,23 @@ void MultiTransport::getMixDiffCoeffs(doublereal* const d)
}
}
void MultiTransport::updateTransport_T()
void MultiTransport::update_T()
{
if (m_temp == m_thermo->temperature()) {
return;
}
m_temp = m_thermo->temperature();
m_logt = log(m_temp);
m_kbt = Boltzmann * m_temp;
m_sqrt_t = sqrt(m_temp);
m_t14 = sqrt(m_sqrt_t);
m_t32 = m_temp * m_sqrt_t;
m_sqrt_kbt = sqrt(Boltzmann*m_temp);
// compute powers of log(T)
m_polytempvec[0] = 1.0;
m_polytempvec[1] = m_logt;
m_polytempvec[2] = m_logt*m_logt;
m_polytempvec[3] = m_logt*m_logt*m_logt;
m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt;
GasTransport::update_T();
// temperature has changed, so polynomial fits will need to be
// redone, and the L matrix reevaluated.
m_visc_ok = false;
m_spvisc_ok = false;
m_diff_ok = false;
m_abc_ok = false;
m_lmatrix_soln_ok = false;
m_l0000_ok = false;
}
void MultiTransport::updateTransport_C()
void MultiTransport::update_C()
{
// signal that concentration-dependent quantities will need to
// be recomputed before use, and update the local mole
@ -828,7 +723,7 @@ void MultiTransport::updateDiff_T()
if (m_diff_tlast == m_thermo->temperature()) {
return;
}
updateTransport_T();
update_T();
// evaluate binary diffusion coefficients at unit pressure
size_t ic = 0;
@ -854,68 +749,13 @@ void MultiTransport::updateDiff_T()
m_diff_tlast = m_thermo->temperature();
}
void MultiTransport::updateSpeciesViscosities_T()
{
if (m_spvisc_tlast == m_thermo->temperature()) {
return;
}
updateTransport_T();
if (m_mode == CK_Mode) {
for (size_t k = 0; k < m_nsp; k++) {
m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
m_sqvisc[k] = sqrt(m_visc[k]);
}
} else {
for (size_t k = 0; k < m_nsp; k++) {
//m_visc[k] = m_sqrt_t*dot5(m_polytempvec, m_visccoeffs[k]);
// the polynomial fit is done for sqrt(visc/sqrt(T))
m_sqvisc[k] = m_t14*dot5(m_polytempvec, m_visccoeffs[k]);
m_visc[k] = (m_sqvisc[k]*m_sqvisc[k]);
}
}
m_spvisc_ok = true;
m_spvisc_tlast = m_thermo->temperature();
}
void MultiTransport::updateViscosity_T()
{
if (m_visc_tlast == m_thermo->temperature()) {
return;
}
doublereal vratiokj, wratiojk, factor1;
updateSpeciesViscosities_T();
// see Eq. (9-5.15) of Reid, Prausnitz, and Poling
for (size_t j = 0; j < m_nsp; j++) {
for (size_t k = j; k < m_nsp; k++) {
vratiokj = m_visc[k]/m_visc[j];
wratiojk = m_mw[j]/m_mw[k];
//rootwjk = sqrt(wratiojk);
//factor1 = 1.0 + sqrt(vratiokj * rootwjk);
//m_phi(k,j) = factor1*factor1 /
// (SqrtEight * sqrt(1.0 + m_mw[k]/m_mw[j]));
//m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
// Note that m_wratjk(k,j) holds the square root of
// m_wratjk(j,k)!
factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
m_phi(k,j) = factor1*factor1 /
(SqrtEight * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
}
m_visc_ok = true;
m_visc_tlast = m_thermo->temperature();
}
void MultiTransport::updateThermal_T()
{
if (m_thermal_tlast == m_thermo->temperature()) {
return;
}
// we need species viscosities and binary diffusion coefficients
updateSpeciesViscosities_T();
updateSpeciesViscosities();
updateDiff_T();
// evaluate polynomial fits for A*, B*, C*