diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index c37cc00b1..67587806f 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -23,7 +23,7 @@ using namespace std; * Mole fractions below MIN_X will be set to MIN_X when computing * transport properties. */ -#define MIN_X 1.e-20 +#define MIN_X 1.e-14 namespace Cantera { @@ -443,6 +443,35 @@ namespace Cantera { } } + void LiquidTransport::getSpeciesDiffusiveMassFluxes(doublereal* const fluxes) { + int n, k; + + update_temp(); + update_conc(); + + + getMixDiffCoeffs(DATA_PTR(m_spwork)); + + + const array_fp& mw = m_thermo->molecularWeights(); + const doublereal* const y = m_thermo->massFractions(); + const doublereal rhon = m_thermo->molarDensity(); + // Unroll wrt ndim + vector_fp sum(m_nDim,0.0); + for (n = 0; n < m_nDim; n++) { + for (k = 0; k < m_nsp; k++) { + fluxes[n*m_nsp + k] = -rhon * mw[k] * m_spwork[k] * m_Grad_X[n*m_nsp + k]; + sum[n] += fluxes[n*m_nsp + k]; + } + } + // add correction flux to enforce sum to zero + for (n = 0; n < m_nDim; n++) { + for (k = 0; k < m_nsp; k++) { + fluxes[n*m_nsp + k] -= y[k]*sum[n]; + } + } + } + /** * Mixture-averaged diffusion coefficients [m^2/s]. * @@ -463,24 +492,28 @@ namespace Cantera { int k, j; doublereal mmw = m_thermo->meanMolecularWeight(); - doublereal sumxw = 0.0, sum2; - doublereal p = m_press; + doublereal sumxw_tran = 0.0; + doublereal sum2; + if (m_nsp == 1) { - d[0] = m_bdiff(0,0) / p; + d[0] = m_bdiff(0,0); } else { - for (k = 0; k < m_nsp; k++) sumxw += m_molefracs[k] * m_mw[k]; + for (k = 0; k < m_nsp; k++) { + sumxw_tran += m_molefracs_tran[k] * m_mw[k]; + } for (k = 0; k < m_nsp; k++) { sum2 = 0.0; for (j = 0; j < m_nsp; j++) { if (j != k) { - sum2 += m_molefracs[j] / m_bdiff(j,k); + sum2 += m_molefracs_tran[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); - } + // Because we use m_molefracs_tran, sum2 must be positive definate + // if (sum2 <= 0.0) { + // d[k] = m_bdiff(k,k); + // } else { + d[k] = (sumxw_tran - m_molefracs_tran[k] * m_mw[k])/(mmw * sum2); + // } } } } @@ -569,12 +602,15 @@ namespace Cantera { if (iStateNew != m_iStateMF) { qReturn = false; m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); - m_thermo->getMoleFractions(DATA_PTR(m_concentrations)); + m_thermo->getConcentrations(DATA_PTR(m_concentrations)); double ctot = 0.0; for (int k = 0; k < m_nsp; k++) { - m_molefracs[k] = fmaxx(MIN_X, m_molefracs[k]); + m_molefracs[k] = fmaxx(0.0, m_molefracs[k]); + m_molefracs_tran[k] = fmaxx(MIN_X, m_molefracs[k]); ctot += 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]); @@ -597,6 +633,12 @@ namespace Cantera { // We formulate the directional derivative + /* + * We only calculate the change in ac due to composition. + * The pressure and the temperature are taken care of in + * other parts of the expression. + * + */ void LiquidTransport::update_Grad_lnAC() { int k; @@ -632,25 +674,12 @@ namespace Cantera { } for (k = 0; k < m_nsp; k++) { - m_Grad_lnAC[m_nsp * a + k] = sum * (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k])) / mag; + m_Grad_lnAC[m_nsp * a + k] = + sum * (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k])) / mag; } } m_thermo->setMoleFractions(DATA_PTR(m_molefracs)); - double Tbase = m_thermo->temperature(); - double T_new = Tbase - 1.0E-6; - m_thermo->setTemperature(T_new); - - m_thermo->getActivityCoefficients(DATA_PTR(lnActCoeffMolarDelta_)); - double *dlnActCoeffdT = &Xdelta_[0]; - for (k = 0; k < m_nsp; k++) { - dlnActCoeffdT[k] == (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k]))/(-1.0E-6); - } - for (int a = 0; a < m_nDim; a++) { - for (k = 0; k < m_nsp; k++) { - m_Grad_lnAC[m_nsp * a + k] += dlnActCoeffdT[k] * m_Grad_T[a]; - } - } - m_thermo->setTemperature(Tbase); + } /************************************************************************* @@ -780,11 +809,7 @@ namespace Cantera { //! grab a local copy of the molecular weights const vector_fp& M = m_thermo->molecularWeights(); - - //! get the mean molecular weight of the mixture - //double M_mix = m_thermo->meanMolecularWeight(); - - + /* * Update the concentrations in the mixture. */ @@ -792,15 +817,24 @@ namespace Cantera { double T = m_thermo->temperature(); - m_thermo->getEntropy_R(DATA_PTR(entropy_R_specSS_)); - m_thermo->getStandardVolumes(DATA_PTR(volume_specSS_)); + + m_thermo->getStandardVolumes(DATA_PTR(volume_specPM_)); m_thermo->getActivityCoefficients(DATA_PTR(actCoeffMolar_)); /* * Calculate the electrochemical potential gradient. This is the * driving force for relative diffusional transport. * - * Here we calculate c_i * grad (mu_i), p. 297 Newman + * Here we calculate + * + * c_i * (grad (mu_i) + S_i grad T - M_i / dens * grad P + * + * This is Eqn. 13-1 p. 318 Newman. The original equation is from + * Hershfeld, Curtis, and Bird. + * + * S_i is the partial molar entropy of species i. This term will cancel + * out a lot of the grad T terms in grad (mu_i), therefore simplifying + * the expression. * * Ok I think there may be many ways to do this. One way is to do it via basis * functions, at the nodes, as a function of the variables in the problem. @@ -808,15 +842,12 @@ namespace Cantera { * For calculation of molality based thermo systems, we current get * the molar based values. This may change. * - * */ for (i = 0; i < m_nsp; 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] * GasConstant * entropy_R_specSS_[i] * m_Grad_T[a] - + m_concentrations[i] * volume_specSS_[i] * m_Grad_P[a] - + m_concentrations[i] * GasConstant * m_Grad_T[a] * log(actCoeffMolar_[i] * m_molefracs[i]) + + 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]; } @@ -829,7 +860,8 @@ namespace Cantera { double lnmnaught = log(mnaught); for (i = 1; i < m_nsp; i++) { for (a = 0; a < VIM; a++) { - m_ck_Grad_mu[a*m_nsp + i] -= m_concentrations[i] * GasConstant * m_Grad_T[a] * lnmnaught; + m_ck_Grad_mu[a*m_nsp + i] -= + m_concentrations[i] * GasConstant * m_Grad_T[a] * lnmnaught; } } } diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h index ec891205f..d9d90673b 100644 --- a/Cantera/src/transport/LiquidTransport.h +++ b/Cantera/src/transport/LiquidTransport.h @@ -33,6 +33,10 @@ namespace Cantera { const int LVISC_WILKES = 1; const int LVISC_MIXTUREAVG = 2; + const int LDIFF_MIXDIFF_UNCORRECTED = 0; + const int LDIFF_MIXDIFF_FLUXCORRECTED = 1; + const int LDIFF_MULTICOMP_STEFANMAXWELL = 2; + class TransportParams; @@ -302,6 +306,20 @@ namespace Cantera { int ldx, const doublereal* grad_X, int ldf, doublereal* fluxes); + //! Return the species diffusive mass fluxes + /*! + * + * + * + * @param ndim The number of spatial dimensions (1, 2, or 3). + * @param grad_T The temperature gradient (ignored in this model). + * @param ldx Leading dimension of the grad_X array. + * The diffusive mass flux of species \e k is computed from + * + * + */ + virtual void getSpeciesDiffusiveMassFluxes(doublereal* const fluxes); + /** * @param ndim The number of spatial dimensions (1, 2, or 3). * @param grad_T The temperature gradient (ignored in this model). @@ -445,9 +463,13 @@ namespace Cantera { //! Array of Binary Diffusivities /*! - * This has a size equal to nsp x nsp - * It is a symmetric matrix. - * D_ii is undefined. + * Depends on the temperature. We have set the pressure dependence + * to zero for this liquid phase constituitve model + * + * This has a size equal to nsp x nsp + * It is a symmetric matrix. + * D_ii is the self diffusion coefficient. D_ii is not + * needed except for when there is one species in the mixture. * * units m2/sec */ @@ -458,8 +480,8 @@ namespace Cantera { * Viscosity of the species * Length = number of species * - * Depends on the temperature and perhaps pressure, but - * not the species concentrations + * Depends on the temperature. We have set the pressure dependence + * to zero for this liquid phase constituitve model * * controlling update boolean -> m_visc_temp_ok */ @@ -496,11 +518,23 @@ 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. * Update info? * length = m_nsp */ vector_fp m_molefracs; + //! Mole fraction vector + /*! + * 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. + * + * Update info? + * length = m_nsp + */ + vector_fp m_molefracs_tran; + vector_fp Xdelta_; //! Local copy of the concentrations of the species in the phase @@ -512,6 +546,8 @@ namespace Cantera { //! Local copy of the total concentration doublereal concTot_; + doublereal meanMolecularWeight_; + doublereal dens_; //! Local copy of the charge of each species /*! @@ -519,9 +555,8 @@ namespace Cantera { */ vector_fp m_chargeSpecies; - vector_fp entropy_R_specSS_; - - vector_fp volume_specSS_; + + vector_fp volume_specPM_; vector_fp actCoeffMolar_; @@ -601,9 +636,13 @@ namespace Cantera { doublereal m_press; //! Solution of the flux system - Array2D m_flux; + /*! + * This is the mass flux of species k + * in units of kg m-3 s-1. + */ + Array2D m_flux; - //! saved value of the mixture thermal conductivity + //! Saved value of the mixture thermal conductivity doublereal m_lambda; //! Saved value of the mixture viscosity