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

Merged mixture-averaged diffusion coefficient calculations in GasTransport

Browse source
This commit is contained in:
Ray Speth 2012-05-11 15:27:50 +00:00
parent f8308d5853
commit 49b631d346
6 changed files with 156 additions and 318 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -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<vector_fp> 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

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

@ -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<vector_fp> 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<vector_fp> 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;
};

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

@ -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<vector_fp> m_diffcoeffs;
// property values
DenseMatrix m_bdiff;
std::vector<std::vector<int> > m_poly;
std::vector<vector_fp> m_astar_poly;
std::vector<vector_fp> 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;

98
src/transport/GasTransport.cpp
View file

@ -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);
}
}
}
}
}

122
src/transport/MixTransport.cpp
View file

@ -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 {
}
//====================================================================================================================
}

112
src/transport/MultiTransport.cpp
View file

@ -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;