From 49b631d34618bfe40e42807ea5031c78f452a89d Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 11 May 2012 15:27:50 +0000 Subject: [PATCH] Merged mixture-averaged diffusion coefficient calculations in GasTransport --- include/cantera/transport/GasTransport.h | 57 ++++++++++ include/cantera/transport/MixTransport.h | 64 ----------- include/cantera/transport/MultiTransport.h | 21 ---- src/transport/GasTransport.cpp | 98 ++++++++++++++++- src/transport/MixTransport.cpp | 122 --------------------- src/transport/MultiTransport.cpp | 112 +------------------ 6 files changed, 156 insertions(+), 318 deletions(-) diff --git a/include/cantera/transport/GasTransport.h b/include/cantera/transport/GasTransport.h index c9eef47bd..6721a3b2b 100644 --- a/include/cantera/transport/GasTransport.h +++ b/include/cantera/transport/GasTransport.h @@ -47,6 +47,35 @@ public: std::copy(m_visc.begin(), m_visc.end(), visc); } + //! Returns the matrix of binary diffusion coefficients. + /*! + * d[ld*j + i] = rp * m_bdiff(i,j); + * + * @param ld offset of rows in the storage + * @param d output vector of diffusion coefficients. Units of m**2 / s + */ + virtual void getBinaryDiffCoeffs(const size_t ld, doublereal* const d); + + //! Returns the Mixture-averaged diffusion coefficients [m^2/s]. + /*! + * Returns the mixture averaged diffusion coefficients for a gas, + * appropriate for calculating the mass averaged diffusive flux with respect + * to the mass averaged velocity using gradients of the mole fraction. + * Note, for the single species case or the pure fluid case the routine + * returns the self-diffusion coefficient. This is needed to avoid a Nan + * result in the formula below. + * + * This is Eqn. 12.180 from "Chemically Reacting Flow" + * + * \f[ + * D_{km}' = \frac{\left( \bar{M} - X_k M_k \right)}{ \bar{\qquad M \qquad } } {\left( \sum_{j \ne k} \frac{X_j}{D_{kj}} \right) }^{-1} + * \f] + * + * @param[out] d Vector of mixture diffusion coefficients, \f$ D_{km}' \f$ , + * for each species (m^2/s). length m_nsp + */ + virtual void getMixDiffCoeffs(doublereal* const d); + protected: GasTransport(ThermoPhase* thermo=0); @@ -73,6 +102,12 @@ protected: //! of pressure. virtual void updateSpeciesViscosities(); + //! Update the binary diffusion coefficients + /*! + * These are evaluated from the polynomial fits of the temperature at the unit pressure of 1 Pa. + */ + void updateDiff_T(); + //! Vector of species mole fractions. These are processed so that all mole //! fractions are >= MIN_X. Length = m_kk. vector_fp m_molefracs; @@ -89,6 +124,9 @@ protected: //! Update boolean for the species viscosities bool m_spvisc_ok; + //! Update boolean for the binary diffusivities at unit pressure + bool m_bindiff_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; @@ -154,6 +192,25 @@ protected: //! Current value of temperature to the 3/2 power doublereal m_t32; + + //! 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: + * + * int ic = 0; + * for (i = 0; i < m_nsp; i++) { + * for (j = i; j < m_nsp; j++) { + * ic++; + * } + * } + */ + std::vector m_diffcoeffs; + + //! Matrix of binary diffusion coefficients at the reference pressure and the current temperature + //! Size is nsp x nsp. + DenseMatrix m_bdiff; }; } // namespace Cantera diff --git a/include/cantera/transport/MixTransport.h b/include/cantera/transport/MixTransport.h index 0e4445773..72fb910fa 100644 --- a/include/cantera/transport/MixTransport.h +++ b/include/cantera/transport/MixTransport.h @@ -143,37 +143,6 @@ public: */ virtual doublereal thermalConductivity(); - //! Returns the matrix of binary diffusion coefficients. - /*! - * - * d[ld*j + i] = rp * m_bdiff(i,j); - * - * units of m**2 / s - * - * @param ld offset of rows in the storage - * @param d output vector of diffusion coefficients - */ - virtual void getBinaryDiffCoeffs(const size_t ld, doublereal* const d); - - //! Returns the Mixture-averaged diffusion coefficients [m^2/s]. - /*! - * Returns the mixture averaged diffusion coefficients for a gas, appropriate for calculating the - * mass averaged diffusive flux with respect to the mass averaged velocity using gradients of the - * mole fraction. - * Note, for the single species case or the pure fluid case the routine returns the self-diffusion coefficient. - * This is need to avoid a Nan result in the formula below. - * - * This is Eqn. 12.180 from "Chemically Reacting Flow" - * - * \f[ - * D_{km}' = \frac{\left( \bar{M} - X_k M_k \right)}{ \bar{\qquad M \qquad } } {\left( \sum_{j \ne k} \frac{X_j}{D_{kj}} \right) }^{-1} - * \f] - * - * @param d Output Vector of mixture diffusion coefficients, \f$ D_{km}' \f$ , for each species (m^2/s). - * length m_nsp - */ - virtual void getMixDiffCoeffs(doublereal* const d); - //! Get the Electrical mobilities (m^2/V/s). /*! * This function returns the mobilities. In some formulations @@ -271,13 +240,6 @@ private: */ void updateCond_T(); - //! Update the binary diffusion coefficients - /*! - * These are evaluated from the polynomial fits of the temperature at the unit pressure of 1 Pa. - */ - void updateDiff_T(); - - // --------- Member Data ------------- private: @@ -288,28 +250,6 @@ private: */ std::vector m_condcoeffs; - //! 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: - * - * int ic = 0; - * for (i = 0; i < m_nsp; i++) { - * for (j = i; j < m_nsp; j++) { - * ic++; - * } - * } - * - */ - std::vector m_diffcoeffs; - - //! Matrix of binary diffusion coefficients at the reference pressure and the current temperature - /*! - * Size is nsp x nsp - */ - DenseMatrix m_bdiff; - //! vector of species thermal conductivities (W/m /K) /*! * These are used in wilke's rule to calculate the viscosity of the solution @@ -324,9 +264,6 @@ private: */ doublereal m_lambda; - //! Update boolean for the binary diffusivities at unit pressure - bool m_bindiff_ok; - //! Update boolean for the species thermal conductivities bool m_spcond_ok; @@ -395,7 +332,6 @@ private: */ vector_fp m_zrot; - //! Debug flag - turns on more printing bool m_debug; }; diff --git a/include/cantera/transport/MultiTransport.h b/include/cantera/transport/MultiTransport.h index f9f10bb15..6aed8c131 100644 --- a/include/cantera/transport/MultiTransport.h +++ b/include/cantera/transport/MultiTransport.h @@ -116,18 +116,8 @@ public: virtual doublereal thermalConductivity(); - virtual void getBinaryDiffCoeffs(const size_t ld, doublereal* const d); virtual void getMultiDiffCoeffs(const size_t ld, doublereal* const d); - //! Although this class implements a multicomponent diffusion - //! model, it is convenient to be able to compute - //! mixture-averaged diffusion coefficients too. - /*! - * @param d Mixture averaged diffusion coefficients - * Length = m_msp, units = m2/sec - */ - virtual void getMixDiffCoeffs(doublereal* const d); - //! Get the species diffusive mass fluxes wrt to the mass averaged velocity, //! given the gradients in mole fraction and temperature /*! @@ -230,24 +220,14 @@ protected: //! conductivity and thermal diffusion coefficients. void updateThermal_T(); - //! Update the binary diffusion coefficients. - //! These are evaluated from the polynomial fits at unit pressure (1 Pa). - void updateDiff_T(); - private: - doublereal m_diff_tlast; doublereal m_thermal_tlast; doublereal m_tmin; doublereal m_tmax; - // polynomial fits - std::vector m_diffcoeffs; - // property values - DenseMatrix m_bdiff; - std::vector > m_poly; std::vector m_astar_poly; std::vector m_bstar_poly; @@ -297,7 +277,6 @@ private: void correctBinDiffCoeffs(); //! Boolean indicating viscosity is up to date - bool m_diff_ok; bool m_abc_ok; bool m_l0000_ok; bool m_lmatrix_soln_ok; diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp index d470b92d3..c1d6dee14 100644 --- a/src/transport/GasTransport.cpp +++ b/src/transport/GasTransport.cpp @@ -10,6 +10,7 @@ GasTransport::GasTransport(ThermoPhase* thermo) : m_visc_ok(false), m_viscwt_ok(false), m_spvisc_ok(false), + m_bindiff_ok(false), m_mode(0), m_phi(0,0), m_spwork(0), @@ -26,7 +27,9 @@ GasTransport::GasTransport(ThermoPhase* thermo) : m_sqrt_t(0.0), m_logt(0.0), m_t14(0.0), - m_t32(0.0) + m_t32(0.0), + m_diffcoeffs(0), + m_bdiff(0, 0) { } @@ -36,6 +39,7 @@ GasTransport::GasTransport(const GasTransport& right) : m_visc_ok(false), m_viscwt_ok(false), m_spvisc_ok(false), + m_bindiff_ok(false), m_mode(0), m_phi(0,0), m_spwork(0), @@ -52,7 +56,9 @@ GasTransport::GasTransport(const GasTransport& right) : m_sqrt_t(0.0), m_logt(0.0), m_t14(0.0), - m_t32(0.0) + m_t32(0.0), + m_diffcoeffs(0), + m_bdiff(0, 0) { } @@ -63,6 +69,7 @@ GasTransport& GasTransport::operator=(const GasTransport& right) m_visc_ok = right.m_visc_ok; m_viscwt_ok = right.m_viscwt_ok; m_spvisc_ok = right.m_spvisc_ok; + m_bindiff_ok = right.m_bindiff_ok; m_mode = right.m_mode; m_phi = right.m_phi; m_spwork = right.m_spwork; @@ -79,6 +86,8 @@ GasTransport& GasTransport::operator=(const GasTransport& right) m_logt = right.m_logt; m_t14 = right.m_t14; m_t32 = right.m_t32; + m_diffcoeffs = right.m_diffcoeffs; + m_bdiff = right.m_bdiff; return *this; } @@ -89,10 +98,16 @@ bool GasTransport::initGas(GasTransportParams& tr) 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); @@ -115,6 +130,7 @@ bool GasTransport::initGas(GasTransportParams& tr) m_visc_ok = false; m_viscwt_ok = false; m_spvisc_ok = false; + m_bindiff_ok = false; } void GasTransport::update_T(void) { @@ -137,6 +153,7 @@ void GasTransport::update_T(void) { m_visc_ok = false; m_spvisc_ok = false; m_viscwt_ok = false; + m_bindiff_ok = false; } doublereal GasTransport::viscosity() @@ -203,4 +220,81 @@ void GasTransport::updateSpeciesViscosities() m_spvisc_ok = true; } +void GasTransport::updateDiff_T() +{ + // evaluate binary diffusion coefficients at unit pressure + size_t ic = 0; + if (m_mode == CK_Mode) { + for (size_t i = 0; i < m_nsp; i++) { + for (size_t j = i; j < m_nsp; j++) { + m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic])); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } else { + for (size_t i = 0; i < m_nsp; i++) { + for (size_t j = i; j < m_nsp; j++) { + m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec, + m_diffcoeffs[ic]); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } + m_bindiff_ok = true; +} + +void GasTransport::getBinaryDiffCoeffs(const size_t ld, doublereal* const d) +{ + update_T(); + // if necessary, evaluate the binary diffusion coefficients from the polynomial fits + if (!m_bindiff_ok) { + updateDiff_T(); + } + if (ld < m_nsp) { + throw CanteraError(" MixTransport::getBinaryDiffCoeffs()", "ld is too small"); + } + doublereal rp = 1.0/m_thermo->pressure(); + for (size_t i = 0; i < m_nsp; i++) + for (size_t j = 0; j < m_nsp; j++) { + d[ld*j + i] = rp * m_bdiff(i,j); + } +} + +void GasTransport::getMixDiffCoeffs(doublereal* const d) +{ + update_T(); + update_C(); + + // update the binary diffusion coefficients if necessary + if (!m_bindiff_ok) { + updateDiff_T(); + } + + doublereal mmw = m_thermo->meanMolecularWeight(); + doublereal sumxw = 0.0, sum2; + doublereal p = m_thermo->pressure(); + if (m_nsp == 1) { + d[0] = m_bdiff(0,0) / p; + } else { + for (size_t k = 0; k < m_nsp; k++) { + sumxw += m_molefracs[k] * m_mw[k]; + } + for (size_t k = 0; k < m_nsp; k++) { + sum2 = 0.0; + for (size_t j = 0; j < m_nsp; j++) { + if (j != k) { + sum2 += m_molefracs[j] / m_bdiff(j,k); + } + } + if (sum2 <= 0.0) { + d[k] = m_bdiff(k,k) / p; + } else { + d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2); + } + } + } +} + } diff --git a/src/transport/MixTransport.cpp b/src/transport/MixTransport.cpp index eee1f2866..effe353d7 100644 --- a/src/transport/MixTransport.cpp +++ b/src/transport/MixTransport.cpp @@ -26,15 +26,11 @@ using namespace std; namespace Cantera { - //==================================================================================================================== MixTransport::MixTransport() : m_condcoeffs(0), - m_diffcoeffs(0), - m_bdiff(0, 0), m_cond(0), m_lambda(0.0), - m_bindiff_ok(false), m_spcond_ok(false), m_condmix_ok(false), m_eps(0), @@ -50,11 +46,8 @@ MixTransport::MixTransport() : MixTransport::MixTransport(const MixTransport& right) : GasTransport(right), m_condcoeffs(0), - m_diffcoeffs(0), - m_bdiff(0, 0), m_cond(0), m_lambda(0.0), - m_bindiff_ok(false), m_spcond_ok(false), m_condmix_ok(false), m_eps(0), @@ -83,11 +76,8 @@ MixTransport& MixTransport::operator=(const MixTransport& right) GasTransport::operator=(right); m_condcoeffs = right.m_condcoeffs; - m_diffcoeffs = right.m_diffcoeffs; - m_bdiff = right.m_bdiff; m_cond = right.m_cond; m_lambda = right.m_lambda; - m_bindiff_ok = right.m_bindiff_ok; m_spcond_ok = right.m_spcond_ok; m_condmix_ok = right.m_condmix_ok; m_eps = right.m_eps; @@ -122,13 +112,10 @@ bool MixTransport::initGas(GasTransportParams& tr) GasTransport::initGas(tr); // copy polynomials and parameters into local storage - m_visccoeffs = tr.visccoeffs; m_condcoeffs = tr.condcoeffs; - m_diffcoeffs = tr.diffcoeffs; m_zrot = tr.zrot; m_crot = tr.crot; - m_mode = tr.mode_; m_diam = tr.diam; m_eps = tr.eps; m_alpha = tr.alpha; @@ -138,7 +125,6 @@ bool MixTransport::initGas(GasTransportParams& tr) } m_cond.resize(m_nsp); - m_bdiff.resize(m_nsp, m_nsp); // set flags all false m_spcond_ok = false; @@ -147,33 +133,6 @@ bool MixTransport::initGas(GasTransportParams& tr) return true; } -//==================================================================================================================== -// Returns the matrix of binary diffusion coefficients. -/* - * - * d[ld*j + i] = rp * m_bdiff(i,j); - * - * units of m**2 / s - * - * @param ld offset of rows in the storage - * @param d output vector of diffusion coefficients - */ -void MixTransport::getBinaryDiffCoeffs(const size_t ld, doublereal* const d) -{ - update_T(); - // if necessary, evaluate the binary diffusion coefficients from the polynomial fits - if (!m_bindiff_ok) { - updateDiff_T(); - } - if (ld < m_nsp) { - throw CanteraError(" MixTransport::getBinaryDiffCoeffs()", "ld is too small"); - } - doublereal rp = 1.0/pressure_ig(); - for (size_t i = 0; i < m_nsp; i++) - for (size_t j = 0; j < m_nsp; j++) { - d[ld*j + i] = rp * m_bdiff(i,j); - } -} //=================================================================================================================== void MixTransport::getMobilities(doublereal* const mobil) { @@ -290,51 +249,6 @@ void MixTransport::getSpeciesFluxes(size_t ndim, } } } -//=========================================================================================================== -// Mixture-averaged diffusion coefficients [m^2/s]. -/* - * Returns the mixture averaged diffusion coefficients for a gas. - * Note, for the single species case or the pure fluid case the routine returns the self-diffusion coefficient. - * This is need to avoid a Nan result in the formula - * below. - * - * @param d Output Vector of diffusion coefficients for each species (m^2/s) - * length m_nsp - */ -void MixTransport::getMixDiffCoeffs(doublereal* const d) -{ - update_T(); - update_C(); - - // update the binary diffusion coefficients if necessary - if (!m_bindiff_ok) { - updateDiff_T(); - } - - doublereal mmw = m_thermo->meanMolecularWeight(); - doublereal sumxw = 0.0, sum2; - doublereal p = pressure_ig(); - if (m_nsp == 1) { - d[0] = m_bdiff(0,0) / p; - } else { - for (size_t k = 0; k < m_nsp; k++) { - sumxw += m_molefracs[k] * m_mw[k]; - } - for (size_t k = 0; k < m_nsp; k++) { - sum2 = 0.0; - for (size_t j = 0; j < m_nsp; j++) { - if (j != k) { - sum2 += m_molefracs[j] / m_bdiff(j,k); - } - } - if (sum2 <= 0.0) { - d[k] = m_bdiff(k,k) / p; - } else { - d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2); - } - } - } -} //=========================================================================================================== /* @@ -399,41 +313,6 @@ void MixTransport::updateCond_T() m_spcond_ok = true; m_condmix_ok = false; } -//==================================================================================================================== -/* - * Update the binary diffusion coefficients. These are evaluated - * from the polynomial fits at unit pressure (1 Pa). - */ -void MixTransport::updateDiff_T() -{ - - // evaluate binary diffusion coefficients at unit pressure - size_t ic = 0; - if (m_mode == CK_Mode) { - for (size_t i = 0; i < m_nsp; i++) { - for (size_t j = i; j < m_nsp; j++) { - m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic])); - m_bdiff(j,i) = m_bdiff(i,j); - ic++; - } - } - } else { - for (size_t i = 0; i < m_nsp; i++) { - for (size_t j = i; j < m_nsp; j++) { - m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec, - m_diffcoeffs[ic]); - m_bdiff(j,i) = m_bdiff(i,j); - ic++; - } - } - } - m_bindiff_ok = true; -} -//==================================================================================================================== -/* - * Update the pure-species viscosities. - */ - //==================================================================================================================== /* @@ -461,4 +340,3 @@ struct GasTransportData MixTransport::getGasTransportData(int kSpecies) const { } //==================================================================================================================== } - diff --git a/src/transport/MultiTransport.cpp b/src/transport/MultiTransport.cpp index 0c7f1faaf..6e5987dbc 100644 --- a/src/transport/MultiTransport.cpp +++ b/src/transport/MultiTransport.cpp @@ -94,8 +94,6 @@ bool MultiTransport::initGas(GasTransportParams& tr) // copy polynomials and parameters into local storage m_poly = tr.poly; - m_visccoeffs = tr.visccoeffs; - m_diffcoeffs = tr.diffcoeffs; m_astar_poly = tr.astar_poly; m_bstar_poly = tr.bstar_poly; m_cstar_poly = tr.cstar_poly; @@ -103,7 +101,6 @@ bool MultiTransport::initGas(GasTransportParams& tr) m_zrot = tr.zrot; m_crot = tr.crot; m_epsilon = tr.epsilon; - m_mode = tr.mode_; m_diam = tr.diam; m_eps = tr.eps; m_alpha = tr.alpha; @@ -123,21 +120,16 @@ bool MultiTransport::initGas(GasTransportParams& tr) m_cinternal.resize(m_nsp); - m_bdiff.resize(m_nsp, m_nsp); - - //m_poly.resize(m_nsp); m_om22.resize(m_nsp, m_nsp); m_astar.resize(m_nsp, m_nsp); m_bstar.resize(m_nsp, m_nsp); m_cstar.resize(m_nsp, m_nsp); // set flags all false - m_diff_ok = false; m_abc_ok = false; m_l0000_ok = false; m_lmatrix_soln_ok = false; - m_diff_tlast = 0.0; m_thermal_tlast = 0.0; // use LU decomposition by default @@ -178,23 +170,6 @@ bool MultiTransport::initGas(GasTransportParams& tr) //==================================================================================================================== -/******************* binary diffusion coefficients **************/ - -void MultiTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d) -{ - // if necessary, evaluate the binary diffusion coefficients - // from the polynomial fits - updateDiff_T(); - - doublereal p = pressure_ig(); - doublereal rp = 1.0/p; - for (size_t i = 0; i < m_nsp; i++) - for (size_t j = 0; j < m_nsp; j++) { - d[ld*j + i] = rp * m_bdiff(i,j); - } -} - - /****************** thermal conductivity **********************/ /** @@ -231,10 +206,7 @@ void MultiTransport::getThermalDiffCoeffs(doublereal* const dt) */ void MultiTransport::solveLMatrixEquation() { - - // if T has changed, update the temperature-dependent - // properties. - + // if T has changed, update the temperature-dependent properties. updateThermal_T(); update_C(); @@ -271,7 +243,6 @@ void MultiTransport::solveLMatrixEquation() } // evaluate the submatrices of the L matrix - m_Lmatrix.resize(3*m_nsp, 3*m_nsp, 0.0); eval_L0000(DATA_PTR(m_molefracs)); @@ -284,7 +255,6 @@ void MultiTransport::solveLMatrixEquation() eval_L0110(); eval_L0101(DATA_PTR(m_molefracs)); - // Solve it using GMRES or LU decomposition. The last solution // in m_a should provide a good starting guess, so convergence // should be fast. @@ -336,15 +306,11 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ const doublereal* const grad_X, int ldf, doublereal* const fluxes) { - // update the binary diffusion coefficients if necessary updateDiff_T(); - doublereal sum; - // If any component of grad_T is non-zero, then get the // thermal diffusion coefficients - bool addThermalDiffusion = false; for (size_t i = 0; i < ndim; i++) { if (grad_T[i] != 0.0) { @@ -359,7 +325,7 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ doublereal rho = m_thermo->density(); for (size_t i = 0; i < m_nsp; i++) { - sum = 0.0; + double sum = 0.0; for (size_t j = 0; j < m_nsp; j++) { m_aa(i,j) = m_molefracs[j]*m_molefracs[i]/m_bdiff(i,j); sum += m_aa(i,j); @@ -381,7 +347,6 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ // set the matrix elements in this row to the mass fractions, // and set the entry in gradx to zero - for (size_t j = 0; j < m_nsp; j++) { m_aa(jmax,j) = y[j]; } @@ -424,7 +389,6 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ throw CanteraError("MultiTransport::getSpeciesFluxes", "Error in DGETRS"); - size_t offset; doublereal pp = pressure_ig(); @@ -468,7 +432,6 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* state2, doublereal delta, doublereal* fluxes) { - double* x1 = DATA_PTR(m_spwork1); double* x2 = DATA_PTR(m_spwork2); double* x3 = DATA_PTR(m_spwork3); @@ -483,7 +446,6 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s double t2 = state2[0]; m_thermo->getMoleFractions(x2); - // double p = 0.5*(p1 + p2); double t = 0.5*(state1[0] + state2[0]); @@ -555,7 +517,6 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s throw CanteraError("MultiTransport::getMassFluxes", "Error in DGETRF. Info = "+int2str(info)); - doublereal pp = pressure_ig(); // multiply diffusion velocities by rho * Y_k to create @@ -646,38 +607,6 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d) } //==================================================================================================================== -void MultiTransport::getMixDiffCoeffs(doublereal* const d) -{ - // update the mole fractions - update_C(); - - // update the binary diffusion coefficients if necessary - updateDiff_T(); - - doublereal mmw = m_thermo->meanMolecularWeight(); - doublereal sumxw = 0.0, sum2; - doublereal p = pressure_ig(); - if (m_nsp == 1) { - d[0] = m_bdiff(0,0) / p; - } else { - for (size_t k = 0; k < m_nsp; k++) { - sumxw += m_molefracs[k] * m_mw[k]; - } - for (size_t k = 0; k < m_nsp; k++) { - sum2 = 0.0; - for (size_t j = 0; j < m_nsp; j++) { - if (j != k) { - sum2 += m_molefracs[j] / m_bdiff(j,k); - } - } - if (sum2 <= 0.0) { - d[k] = m_bdiff(k,k) / p; - } else { - d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2); - } - } - } -} void MultiTransport::update_T() { @@ -689,7 +618,6 @@ void MultiTransport::update_T() // temperature has changed, so polynomial fits will need to be // redone, and the L matrix reevaluated. - m_diff_ok = false; m_abc_ok = false; m_lmatrix_soln_ok = false; m_l0000_ok = false; @@ -704,7 +632,6 @@ void MultiTransport::update_C() m_lmatrix_soln_ok = false; m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); - // add an offset to avoid a pure species condition // (check - this may be unnecessary) for (size_t k = 0; k < m_nsp; k++) { @@ -718,37 +645,6 @@ void MultiTransport::update_C() * *************************************************************************/ -void MultiTransport::updateDiff_T() -{ - if (m_diff_tlast == m_thermo->temperature()) { - return; - } - update_T(); - - // evaluate binary diffusion coefficients at unit pressure - size_t ic = 0; - if (m_mode == CK_Mode) { - for (size_t i = 0; i < m_nsp; i++) { - for (size_t j = i; j < m_nsp; j++) { - m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic])); - m_bdiff(j,i) = m_bdiff(i,j); - ic++; - } - } - } else { - for (size_t i = 0; i < m_nsp; i++) { - for (size_t j = i; j < m_nsp; j++) { - m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec, - m_diffcoeffs[ic]); - m_bdiff(j,i) = m_bdiff(i,j); - ic++; - } - } - } - m_diff_ok = true; - m_diff_tlast = m_thermo->temperature(); -} - void MultiTransport::updateThermal_T() { if (m_thermal_tlast == m_thermo->temperature()) { @@ -784,9 +680,7 @@ void MultiTransport::updateThermal_T() } m_abc_ok = true; - // evaluate the temperature-dependent rotational relaxation - // rate - + // evaluate the temperature-dependent rotational relaxation rate doublereal tr, sqtr; for (size_t k = 0; k < m_nsp; k++) { tr = m_eps[k]/ m_kbt;