From f8308d5853e95e4e1516db164278d3d07ec34bcb Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 11 May 2012 15:27:44 +0000 Subject: [PATCH] Merged viscosity calculations from MixTransport and MultiTransport --- include/cantera/transport/GasTransport.h | 161 +++++++++++++++ include/cantera/transport/MixTransport.h | 178 +--------------- include/cantera/transport/MultiTransport.h | 57 +---- src/transport/GasTransport.cpp | 206 ++++++++++++++++++ src/transport/MixTransport.cpp | 230 +-------------------- src/transport/MultiTransport.cpp | 180 +--------------- 6 files changed, 394 insertions(+), 618 deletions(-) create mode 100644 include/cantera/transport/GasTransport.h create mode 100644 src/transport/GasTransport.cpp diff --git a/include/cantera/transport/GasTransport.h b/include/cantera/transport/GasTransport.h new file mode 100644 index 000000000..c9eef47bd --- /dev/null +++ b/include/cantera/transport/GasTransport.h @@ -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 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 diff --git a/include/cantera/transport/MixTransport.h b/include/cantera/transport/MixTransport.h index b2dc7c742..0e4445773 100644 --- a/include/cantera/transport/MixTransport.h +++ b/include/cantera/transport/MixTransport.h @@ -18,16 +18,14 @@ #include // 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 m_visccoeffs; //! Polynomial fits to the thermal conductivity of each species /*! @@ -383,32 +304,12 @@ private: */ std::vector 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 diff --git a/include/cantera/transport/MultiTransport.h b/include/cantera/transport/MultiTransport.h index f2d1f5fbf..f9f10bb15 100644 --- a/include/cantera/transport/MultiTransport.h +++ b/include/cantera/transport/MultiTransport.h @@ -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 m_visccoeffs; - std::vector m_diffcoeffs; - vector_fp m_polytempvec; + std::vector m_diffcoeffs; // property values DenseMatrix m_bdiff; - vector_fp m_visc; - vector_fp m_sqvisc; - - vector_fp m_molefracs; - std::vector > m_poly; std::vector 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(); diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp new file mode 100644 index 000000000..d470b92d3 --- /dev/null +++ b/src/transport/GasTransport.cpp @@ -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; +} + +} diff --git a/src/transport/MixTransport.cpp b/src/transport/MixTransport.cpp index ec6f3c2b2..eee1f2866 100644 --- a/src/transport/MixTransport.cpp +++ b/src/transport/MixTransport.cpp @@ -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(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. diff --git a/src/transport/MultiTransport.cpp b/src/transport/MultiTransport.cpp index b7922f227..0c7f1faaf 100644 --- a/src/transport/MultiTransport.cpp +++ b/src/transport/MultiTransport.cpp @@ -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(*this)); - // m_update_transport_C = m_thermo->installUpdater_C( - // new UpdateTransport_C(*this)); - // m_update_spvisc_T = m_thermo->installUpdater_T( - // new UpdateSpeciesVisc(*this)); - // m_update_visc_T = m_thermo->installUpdater_T( - // new UpdateVisc_T(*this)); - // m_update_diff_T = m_thermo->installUpdater_T( - // new UpdateDiff_T(*this)); - // m_update_thermal_T = m_thermo->installUpdater_T( - // new UpdateThermal_T(*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*