diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index 67587806f..a1e9cc6a5 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -603,18 +603,17 @@ namespace Cantera { qReturn = false; m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); m_thermo->getConcentrations(DATA_PTR(m_concentrations)); - double ctot = 0.0; + concTot_ = 0.0; + concTot_tran_ = 0.0; for (int k = 0; k < m_nsp; k++) { m_molefracs[k] = fmaxx(0.0, m_molefracs[k]); m_molefracs_tran[k] = fmaxx(MIN_X, m_molefracs[k]); - ctot += m_concentrations[k]; + concTot_tran_ += m_molefracs_tran[k]; + concTot_ += m_concentrations[k]; } dens_ = m_thermo->density(); meanMolecularWeight_ = m_thermo->meanMolecularWeight(); - double ctotmin = 0.0; - for (int k = 0; k < m_nsp; k++) { - m_concentrations[k]= fmaxx(ctotmin, m_concentrations[k]); - } + concTot_tran_ *= concTot_; } if (qReturn) { return; @@ -803,7 +802,7 @@ namespace Cantera { */ void LiquidTransport::stefan_maxwell_solve() { int i, j, a; - + doublereal tmp; int VIM = m_nDim; m_B.resize(m_nsp, VIM); //! grab a local copy of the molecular weights @@ -842,18 +841,22 @@ namespace Cantera { * For calculation of molality based thermo systems, we current get * the molar based values. This may change. * + * Note, we have broken the symmetry of the matrix here, due to + * consideratins involving species concentrations going to zero. + * */ for (i = 0; i < m_nsp; i++) { + double xi_denom = m_molefracs_tran[i]; for (a = 0; a < VIM; a++) { m_ck_Grad_mu[a*m_nsp + i] = - m_chargeSpecies[i] * m_concentrations[i] * Faraday * m_Grad_V[a] - + m_concentrations[i] * (volume_specPM_[i] - M[i]/dens_) * m_Grad_P[a] - + m_concentrations[i] * GasConstant * T * m_Grad_lnAC[a*m_nsp+i] / actCoeffMolar_[i] - + concTot_ * GasConstant * T * m_Grad_X[a*m_nsp+i]; + m_chargeSpecies[i] * concTot_ * Faraday * m_Grad_V[a] + + concTot_ * (volume_specPM_[i] - M[i]/dens_) * m_Grad_P[a] + + concTot_ * GasConstant * T * m_Grad_lnAC[a*m_nsp+i] / actCoeffMolar_[i] + + concTot_ * GasConstant * T * m_Grad_X[a*m_nsp+i] / xi_denom; } } - if (m_thermo->activityConvention() == cAC_CONVENTION_MOLALITY ) { + if (m_thermo->activityConvention() == cAC_CONVENTION_MOLALITY) { int iSolvent = 0; double mwSolvent = m_thermo->molecularWeight(iSolvent); double mnaught = mwSolvent/ 1000.; @@ -870,7 +873,7 @@ namespace Cantera { * Just for Note, m_A(i,j) refers to the ith row and jth column. * They are still fortran ordered, so that i varies fastest. */ - switch ( VIM ) { + switch (VIM) { case 1: /* 1-D approximation */ m_B(0,0) = 0.0; for (j = 0; j < m_nsp; j++) { @@ -878,14 +881,12 @@ namespace Cantera { } for (i = 1; i < m_nsp; i++){ m_B(i,0) = m_ck_Grad_mu[i] / (GasConstant * T); + m_A(i,i) = 0.0; for (j = 0; j < m_nsp; j++){ if (j != i) { - m_A(i,j) = m_concentrations[i] * m_concentrations[j]/ - (concTot_ * m_DiffCoeff_StefMax(i,j)); - m_A(j,i) = -m_A(i,j); - } - else if (j == i) { - m_A(i,i) = 0.0; + tmp = m_concentrations[j]/ m_DiffCoeff_StefMax(i,j); + m_A(i,i) += tmp; + m_A(i,j) = - tmp; } } } @@ -901,17 +902,14 @@ namespace Cantera { m_A(0,j) = M[j] * m_concentrations[j]; } for (i = 1; i < m_nsp; i++){ - m_B(i,0) = m_ck_Grad_mu[i] / (GasConstant * T); + m_B(i,0) = m_ck_Grad_mu[i] / (GasConstant * T); m_B(i,1) = m_ck_Grad_mu[m_nsp + i] / (GasConstant * T); - for (j = 0; j < m_nsp; j++){ + m_A(i,i) = 0.0; + for (j = 0; j < m_nsp; j++) { if (j != i) { - m_A(i,j) = m_concentrations[i] * m_concentrations[j]/ - (concTot_ * m_DiffCoeff_StefMax(i,j)); - m_A(j,i) = -m_A(i,j); - - } - else if (j == i) { - m_A(i,i) = 0.0; + tmp = m_concentrations[j] / m_DiffCoeff_StefMax(i,j); + m_A(i,i) += tmp; + m_A(i,j) = - tmp; } } } @@ -930,18 +928,15 @@ namespace Cantera { m_A(0,j) = M[j] * m_concentrations[j]; } for (i = 1; i < m_nsp; i++){ - m_B(i,0) = m_ck_Grad_mu[i] / (GasConstant * T); - m_B(i,1) = m_ck_Grad_mu[m_nsp + i] / (GasConstant * T); + m_B(i,0) = m_ck_Grad_mu[i] / (GasConstant * T); + m_B(i,1) = m_ck_Grad_mu[m_nsp + i] / (GasConstant * T); m_B(i,2) = m_ck_Grad_mu[2*m_nsp + i] / (GasConstant * T); - for (j = 0; j < m_nsp; j++){ + m_A(i,i) = 0.0; + for (j = 0; j < m_nsp; j++) { if (j != i) { - m_A(i,j) = m_concentrations[i] * m_concentrations[j]/ - (concTot_ * m_DiffCoeff_StefMax(i,j)); - m_A(j,i) = -m_A(i,j); - - } - else if (j == i) { - m_A(i,i) = 0.0; + tmp = m_concentrations[j]/ m_DiffCoeff_StefMax(i,j); + m_A(i,i) += tmp; + m_A(i,j) = - tmp; } } } diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h index d9d90673b..e1c72a086 100644 --- a/Cantera/src/transport/LiquidTransport.h +++ b/Cantera/src/transport/LiquidTransport.h @@ -518,17 +518,22 @@ namespace Cantera { //! Local copy of the mole fractions of the species in the phase /*! - * The mole fractions here are assumed to be bounded by 0.0, and 1.0 - * and they are assumed to add up to one. + * The mole fractions here are assumed to be bounded by 0.0 and 1.0 + * and they are assumed to add up to one exactly. This mole + * fraction vector comes from the ThermoPhase object. Derivative + * quantities from this are referred to as bounded. + * * Update info? * length = m_nsp */ vector_fp m_molefracs; - //! Mole fraction vector + //! Non-zero mole fraction vector used in transport property calculations /*! * The mole fractions here are assumed to be bounded by MIN_X and 1.0 - * and they may not be assumed to add up to one. + * and they may not be assumed to add up to one. This + * mole fraction vector is created from the ThermoPhase object. + * Derivative quantities of this use the _tran suffix. * * Update info? * length = m_nsp @@ -539,13 +544,28 @@ namespace Cantera { //! Local copy of the concentrations of the species in the phase /*! + * The concentrations are consistent with the m_molefracs + * vector which is bounded and sums to one. + * * Update info? * length = m_nsp */ vector_fp m_concentrations; - //! Local copy of the total concentration + //! Local copy of the total concentration. + /*! + * This is consistent with the m_concentrations[] and + * m_molefracs[] vector. + */ doublereal concTot_; + + //! Local copy of the total concentration. + /*! + * This is consistent with the x_molefracs_tran vector and + * with the concTot_ number; + */ + doublereal concTot_tran_; + doublereal meanMolecularWeight_; doublereal dens_;