From 433914e4496c9ab02b9894b1364f35a838955714 Mon Sep 17 00:00:00 2001 From: John Hewson Date: Thu, 12 Nov 2009 23:57:47 +0000 Subject: [PATCH] LiquidTransport.cpp and LiquidTransport.h Replaced viscosityModel_ member with m_viscMixModel. Fixed bug in filling m_coeffVisc_Ns, etc., for the Arrhenius model. Failed to use "push_back() to add to vector prior to this. For temperature-dependence option LTR_MODEL_POLY, corrected the polynomial evaluations to correctly account for number of terms in the polynomial. Viscosity mixture evaluations are now up to date with respect to the use of LTR_MIXMODEL_MOLEFRACS and LTR_MIXMODEL_LOG_MOLEFRACS. Added membersto hold activity coefficient info vector_fp actCoeffMolar_; vector_fp lnActCoeffMolarDelta_; Added method getSpeciesFluxesES(...,const doublereal* grad_Phi,doublereal* fluxes)) that takes electrostatic potential gradient in addition to other gradients. This then calls getSpeciesFluxesExt() to do the work. Added new version of method LiquidTransport::update_Grad_lnAC(). The new version calls m_thermo->getdlnActCoeffdlnX( DATA_PTR(grad_lnAC) ); These routines in the thermo object need to be provided. --- Cantera/src/transport/LiquidTransport.cpp | 139 ++++++++++++++++------ Cantera/src/transport/LiquidTransport.h | 62 ++++++++-- 2 files changed, 153 insertions(+), 48 deletions(-) diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index cd196c725..794e6eeda 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -55,7 +55,8 @@ namespace Cantera { m_cond_temp_ok(false), m_cond_mix_ok(false), m_mode(-1000), - m_debug(false) + m_debug(false), + m_nDim(1) { } @@ -81,7 +82,8 @@ namespace Cantera { m_cond_temp_ok(false), m_cond_mix_ok(false), m_mode(-1000), - m_debug(false) + m_debug(false), + m_nDim(1) { /* * Use the assignment operator to do the brunt @@ -126,9 +128,10 @@ namespace Cantera { m_molefracs = right.m_molefracs; m_molefracs_tran = right.m_molefracs_tran; m_concentrations = right.m_concentrations; + m_actCoeff = right.m_actCoeff; + m_Grad_lnAC = right.m_Grad_lnAC; m_chargeSpecies = right.m_chargeSpecies; m_DiffCoeff_StefMax = right.m_DiffCoeff_StefMax; - viscosityModel_ = right.viscosityModel_; m_B = right.m_B; m_A = right.m_A; m_temp = right.m_temp; @@ -201,7 +204,7 @@ namespace Cantera { } else if ( m_viscTempDepType_Ns[k] == LTR_MODEL_ARRHENIUS ) { kentry = ltd.viscCoeffs; //for Arrhenius form, also carry the logarithm of the pre-exponential - kentry[3] = log( kentry[0] ); + kentry.push_back( log( kentry[0] ) ); //should be entry [3] } else if ( m_viscTempDepType_Ns[k] == LTR_MODEL_NOTSET ) { //we might be OK with viscosity not being set so @@ -242,7 +245,7 @@ namespace Cantera { } else if ( m_lambdaTempDepType_Ns[k] == LTR_MODEL_ARRHENIUS ) { kentry = ltd.thermalCondCoeffs; //for Arrhenius form, also carry the logarithm of the pre-exponential - kentry[3] = log( kentry[0] ); + kentry.push_back( log( kentry[0] ) );//should be entry [3] } else if ( m_lambdaTempDepType_Ns[k] == LTR_MODEL_NOTSET ) { throw CanteraError("LiquidTransport::initLiquid", @@ -280,7 +283,7 @@ namespace Cantera { } else if ( m_radiusTempDepType_Ns[k] == LTR_MODEL_ARRHENIUS ) { kentry = ltd.hydroRadiusCoeffs; //for Arrhenius form, also carry the logarithm of the pre-exponential - kentry[3] = log( kentry[0] ); + kentry.push_back( log( kentry[0] ) );//should be entry [3] } else if ( m_radiusTempDepType_Ns[k] == LTR_MODEL_NOTSET ) { throw CanteraError("LiquidTransport::initLiquid", @@ -356,7 +359,7 @@ namespace Cantera { /* * Hydrodynamic radius mixing model rules */ - m_radiusMixModel = tr.model_radius; + m_radiusMixModel = tr.model_hydroradius; m_radius_Aij.resize(m_nsp,m_nsp); m_radius_Aij = tr.radius_Aij; @@ -373,6 +376,8 @@ namespace Cantera { m_molefracs.resize(m_nsp); m_molefracs_tran.resize(m_nsp); m_concentrations.resize(m_nsp); + m_actCoeff.resize(m_nsp); + m_Grad_lnAC.resize(m_nsp); m_spwork.resize(m_nsp); // resize the internal gradient variables @@ -433,16 +438,19 @@ namespace Cantera { /* We still need to implement interaction parameters */ /* This constant viscosity model has no input */ - if (viscosityModel_ == LVISC_CONSTANT) { + if (m_viscMixModel == LTR_MIXMODEL_NOTSET) { - err("constant viscosity not implemented for LiquidTransport."); + err("A viscosity mixing model must be implemented for LiquidTransport."); //return m_viscmix; - } else if (viscosityModel_ == LVISC_AVG_ENERGIES) { + } else if (m_viscMixModel == LTR_MIXMODEL_MOLEFRACS) { - m_viscmix = exp( dot_product(m_logViscSpecies, m_molefracs) ); + m_viscmix = dot_product(m_viscSpecies, m_molefracs) ; + for ( int i = 0; i < m_nsp; i++ ) + for ( int j = 0; j < i; j++ ) + m_viscmix += m_molefracs[i] * m_molefracs[j] * m_visc_Sij(i,j) ; - } else if (viscosityModel_ == LVISC_INTERACTION) { + } else if (m_viscMixModel == LTR_MIXMODEL_LOG_MOLEFRACS) { // log_visc_mix = sum_i (X_i log_visc_i) + sum_i sum_j X_i X_j G_ij double interaction = dot_product(m_logViscSpecies, m_molefracs); @@ -565,21 +573,26 @@ namespace Cantera { m_Grad_V[a] = grad_V[a]; } } - //================================================================================================ + //============================================================== void LiquidTransport::set_Grad_T(const doublereal* const grad_T) { for (int a = 0; a < m_nDim; a++) { m_Grad_T[a] = grad_T[a]; } } - //================================================================================================ + //============================================================== + void LiquidTransport::set_Grad_V(const doublereal* const grad_Phi) { + for (int a = 0; a < m_nDim; a++) { + m_Grad_V[a] = grad_Phi[a]; + } + } + //============================================================== void LiquidTransport::set_Grad_X(const doublereal* const grad_X) { int itop = m_nDim * m_nsp; for (int i = 0; i < itop; i++) { m_Grad_X[i] = grad_X[i]; } - update_Grad_lnAC(); } - //================================================================================================ + //============================================================== /****************** thermal conductivity **********************/ /* * The thermal conductivity is computed from the following mixture rule: @@ -647,6 +660,29 @@ namespace Cantera { getSpeciesFluxesExt(ldf, 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 + * + * \f[ + * \vec{j}_k = -n M_k D_k \nabla X_k. + * \f] + */ + void LiquidTransport::getSpeciesFluxesES(int ndim, + const doublereal* grad_T, + int ldx, + const doublereal* grad_X, + int ldf, + const doublereal* grad_Phi, + doublereal* fluxes) { + set_Grad_T(grad_T); + set_Grad_X(grad_X); + set_Grad_Phi(grad_Phi); + getSpeciesFluxesExt(ldf, fluxes); + } + /** * @param ndim The number of spatial dimensions (1, 2, or 3). * @param grad_T The temperature gradient (ignored in this model). @@ -663,6 +699,7 @@ namespace Cantera { update_T(); update_C(); + update_Grad_lnAC(); getMixDiffCoeffs(DATA_PTR(m_spwork)); @@ -843,10 +880,10 @@ namespace Cantera { * other parts of the expression. * */ + /* void LiquidTransport::update_Grad_lnAC() { int k; - for (int a = 0; a < m_nDim; a++) { // We form the directional derivative double * ma_Grad_X = &m_Grad_X[a*m_nsp]; @@ -885,6 +922,36 @@ namespace Cantera { m_thermo->setMoleFractions(DATA_PTR(m_molefracs)); } + */ + + //! Evaluate the gradient of the activity coefficients + //! as they alter the diffusion coefficient. + /** + * The required quantity is the derivitive of the logarithm of the + * activity coefficient with respect to the derivative of the + * logarithm of the mole fraction (or whatever concentration + * variable we are using to express chemical potential. + * + * Returns the vector over species i: + * \[ + * 1 + \partial \left[ \ln ( \gamma_i ) \right] + * / \partial \left[ \ln ( \X_i ) \right] + * \] + */ + void LiquidTransport::update_Grad_lnAC() { + + int k; + + vector_fp grad_lnAC(m_nsp); + m_thermo->getdlnActCoeffdlnX( DATA_PTR(grad_lnAC) ); + + for (k = 0; k < m_nsp; k++) { + m_Grad_lnAC[k] = 1.0 + grad_lnAC[k]; + std::cout << k << " m_Grad_lnAC = " << m_Grad_lnAC[k] << std::endl; + } + + return; + } /************************************************************************* * @@ -915,12 +982,11 @@ namespace Cantera { - coeffk[2] / m_temp ); } else if ( m_lambdaTempDepType_Ns[k] == LTR_MODEL_POLY ) { - m_lambdaSpecies[k] = coeffk[0] - + coeffk[1] * m_temp - + coeffk[2] * m_temp * m_temp - + coeffk[3] * m_temp * m_temp * m_temp - + coeffk[4] * m_temp * m_temp * m_temp * m_temp; - + double tempN = 1.0; + for ( int i = 0; i < coeffk.size() ; i++ ) { + m_lambdaSpecies[k] += coeffk[i] * tempN; + tempN *= m_temp; + } } else if ( m_lambdaTempDepType_Ns[k] == LTR_MODEL_NOTSET ) { throw CanteraError("LiquidTransport::updateCond_T", "Conductivity Model is not set for species " @@ -1004,12 +1070,11 @@ namespace Cantera { m_viscSpecies[k] = exp( m_logViscSpecies[k] ); } else if ( m_viscTempDepType_Ns[k] == LTR_MODEL_POLY ) { - m_viscSpecies[k] = coeffk[0] - + coeffk[1] * m_temp - + coeffk[2] * m_temp * m_temp - + coeffk[3] * m_temp * m_temp * m_temp - + coeffk[4] * m_temp * m_temp * m_temp * m_temp; - m_logViscSpecies[k] = log( m_viscSpecies[k] ); + double tempN = 1.0; + for ( int i = 0; i < coeffk.size() ; i++ ) { + m_viscSpecies[k] += coeffk[i] * tempN; + tempN *= m_temp; + } } else if ( m_viscTempDepType_Ns[k] == LTR_MODEL_NOTSET ) { throw CanteraError("LiquidTransport::updateViscosity_T", @@ -1058,12 +1123,11 @@ namespace Cantera { - coeffk[2] / m_temp ); } else if ( m_radiusTempDepType_Ns[k] == LTR_MODEL_POLY ) { - m_hydrodynamic_radius[k] = coeffk[0] - + coeffk[1] * m_temp - + coeffk[2] * m_temp * m_temp - + coeffk[3] * m_temp * m_temp * m_temp - + coeffk[4] * m_temp * m_temp * m_temp * m_temp; - + double tempN = 1.0; + for ( int i = 0; i < coeffk.size() ; i++ ) { + m_hydrodynamic_radius[k] += coeffk[i] * tempN; + tempN *= m_temp; + } } else if ( m_radiusTempDepType_Ns[k] == LTR_MODEL_NOTSET ) { throw CanteraError("LiquidTransport::updateHydrodynamicRadius_T", "Hydrodynamic Radius Model is not set for species " @@ -1103,9 +1167,10 @@ namespace Cantera { double T = m_thermo->temperature(); + update_Grad_lnAC() ; m_thermo->getStandardVolumes(DATA_PTR(volume_specPM_)); - m_thermo->getActivityCoefficients(DATA_PTR(actCoeffMolar_)); + m_thermo->getActivityCoefficients(DATA_PTR(m_actCoeff)); /* * Calculate the electrochemical potential gradient. This is the @@ -1138,7 +1203,7 @@ namespace Cantera { m_ck_Grad_mu[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_lnAC[a*m_nsp+i] / m_actCoeff[i] + concTot_ * GasConstant * T * m_Grad_X[a*m_nsp+i] / xi_denom; } } diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h index a8be0e628..7c18037db 100644 --- a/Cantera/src/transport/LiquidTransport.h +++ b/Cantera/src/transport/LiquidTransport.h @@ -314,7 +314,7 @@ namespace Cantera { //! Specify the value of the gradient of the temperature /*! * - * @param grad_V Gradient of the temperature (length num dimensions); + * @param grad_T Gradient of the temperature (length num dimensions); */ virtual void set_Grad_T(const doublereal* const grad_T); @@ -336,7 +336,7 @@ namespace Cantera { * \f[ * \nabla \mu_k = RT \nabla ( \ln X_k ) - * \[ 1 + \nabla ( \ln \gamma_k ) / \nabla ( \ln X_k ) \] + * \left[ 1 + \nabla ( \ln \gamma_k ) / \nabla ( \ln X_k ) \right] * \f] * * The quantity within the square brackets is computed within @@ -369,6 +369,38 @@ namespace Cantera { int ldx, const doublereal* grad_X, int ldf, doublereal* fluxes); + //! Return the species diffusive mass fluxes wrt to + //! the mole averaged velocity, + /** + * @param ndim The number of spatial dimensions (1, 2, or 3). + * @param grad_T The temperature gradient (ignored in this model). + * (length = ndim) + * @param ldx Leading dimension of the grad_X array. + * (usually equal to m_nsp but not always) + * @param grad_X Gradients of the mole fraction + * Flat vector with the m_nsp in the inner loop. + * length = ldx * ndim + * @param ldf Leading dimension of the fluxes array + * (usually equal to m_nsp but not always) + * @param grad_Phi Gradients of the electrostatic potential + * length = ndim + * @param fluxes Output of the diffusive mass fluxes + * Flat vector with the m_nsp in the inner loop. + * length = ldx * ndim + * + * + * The diffusive mass flux of species \e k is computed from + * + * + */ + virtual void getSpeciesFluxesES(int ndim, + const doublereal* grad_T, + int ldx, + const doublereal* grad_X, + int ldf, + const doublereal* grad_Phi, + doublereal* fluxes); + //! Return the species diffusive mass fluxes wrt to //! the mass averaged velocity, /*! @@ -644,6 +676,22 @@ namespace Cantera { */ vector_fp m_Grad_X; + //! Gradient of the logarithm of the activity coefficients + //! with respect to the logarithm of the mole fraction, plus one. + /*! + * This quantity appears in the gradient of the chemical potential. + * It multiplies the gradient of the mole fraction, and in this way + * serves to "modify" the diffusion coefficient. + * + * m_Grad_X[k] = 1 + \partial \left[ \ln ( \gamma_i ) \right] + * / \partial \left[ \ln ( \X_i ) \right] + * + * Note that where "molefraction is used here, whatever + * concentration-related variable applies, so that if + * molality is the concentration variable, the gradient of the + * activity coefficient should be with respect to the molality. + * + */ vector_fp m_Grad_lnAC; //! Internal value of the gradient of the Temperature vector @@ -800,9 +848,7 @@ namespace Cantera { vector_fp volume_specPM_; - vector_fp actCoeffMolar_; - - vector_fp lnActCoeffMolarDelta_; + vector_fp m_actCoeff; //! Stefan-Maxwell Diffusion Coefficients at T, P and C /*! @@ -811,12 +857,6 @@ namespace Cantera { */ DenseMatrix m_DiffCoeff_StefMax; - //! Viscosity model - /*! - * - */ - int viscosityModel_; - //! RHS to the stefan-maxwell equation DenseMatrix m_B;