diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index bf69d190e..aa2f0a3da 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -39,10 +39,6 @@ namespace Cantera { m_iStateMF(-1), m_temp(-1.0), m_logt(0.0), - m_sqrt_t(-1.0), - m_t14(-1.0), - m_t32(-1.0), - m_sqrt_kbt(-1.0), m_press(-1.0), m_lambda(-1.0), m_viscmix(-1.0), @@ -67,10 +63,6 @@ namespace Cantera { m_iStateMF(-1), m_temp(-1.0), m_logt(0.0), - m_sqrt_t(-1.0), - m_t14(-1.0), - m_t32(-1.0), - m_sqrt_kbt(-1.0), m_press(-1.0), m_lambda(-1.0), m_viscmix(-1.0), @@ -101,19 +93,24 @@ namespace Cantera { m_tmax = right.m_tmax; m_mw = right.m_mw; m_visc_A = right.m_visc_A; + m_visc_logA = right.m_visc_logA; m_visc_n = right.m_visc_n; m_visc_Tact = right.m_visc_Tact; + m_visc_Eij = right.m_visc_Eij; + m_visc_Sij = right.m_visc_Sij; m_thermCond_A = right.m_thermCond_A; m_thermCond_n = right.m_thermCond_n; m_thermCond_Tact = right.m_thermCond_Tact; + m_hydrodynamic_radius = right.m_hydrodynamic_radius; m_diffcoeffs = right.m_diffcoeffs; m_Grad_X = right.m_Grad_X; m_Grad_T = right.m_Grad_T; m_Grad_V = right.m_Grad_V; m_ck_Grad_mu = right.m_ck_Grad_mu; m_bdiff = right.m_bdiff; - m_viscSpecies = right.m_viscSpecies; - m_cond = right.m_cond; + m_viscSpecies = right.m_viscSpecies; + m_logViscSpecies = right.m_logViscSpecies; + m_condSpecies = right.m_condSpecies; m_iStateMF = -1; m_molefracs = right.m_molefracs; m_concentrations = right.m_concentrations; @@ -122,13 +119,8 @@ namespace Cantera { viscosityModel_ = right.viscosityModel_; m_B = right.m_B; m_A = right.m_A; - m_eps = right.m_eps; m_temp = right.m_temp; m_logt = right.m_logt; - m_sqrt_t = right.m_sqrt_t; - m_t14 = right.m_t14; - m_t32 = right.m_t32; - m_sqrt_kbt = right.m_sqrt_kbt; m_press = right.m_press; m_flux = right.m_flux; m_lambda = right.m_lambda; @@ -177,23 +169,28 @@ namespace Cantera { m_visc_n = tr.visc_n ; m_visc_Tact = tr.visc_Tact ; + //The following two are not yet filled in LiquidTransportParams + m_visc_Eij = tr.visc_Eij ; + m_visc_Sij = tr.visc_Sij ; + + //save logarithm of pre-exponential for easier computation + m_visc_logA.resize(m_nsp); + for ( int i = 0; i < m_nsp; i++ ) + m_visc_logA[i] = log( m_visc_A[i] ); + m_thermCond_A = tr.thermCond_A ; m_thermCond_n = tr.thermCond_n ; m_thermCond_Tact = tr.thermCond_Tact ; + + m_hydrodynamic_radius = tr.hydroRadius ; + //m_diffcoeffs = tr.diffcoeffs; m_mode = tr.mode_; - m_visc_A.resize(m_nsp); - m_visc_n.resize(m_nsp); - m_visc_Tact.resize(m_nsp); - - m_thermCond_A.resize(m_nsp); - m_thermCond_n.resize(m_nsp); - m_thermCond_Tact.resize(m_nsp); - m_viscSpecies.resize(m_nsp); + m_logViscSpecies.resize(m_nsp); m_condSpecies.resize(m_nsp); m_bdiff.resize(m_nsp, m_nsp); @@ -257,15 +254,24 @@ namespace Cantera { /* We still need to implement interaction parameters */ /* This constant viscosity model has no input */ if (viscosityModel_ == LVISC_CONSTANT) { + err("constant viscosity not implemented for LiquidTransport."); //return m_viscmix; - } else if (viscosityModel_ == LVISC_MIXTUREAVG) { - m_viscmix = dot_product(m_viscSpecies, m_molefracs); + + } else if (viscosityModel_ == LVISC_AVG_ENERGIES) { + + m_viscmix = exp( dot_product(m_logViscSpecies, m_molefracs) ); + } else if (viscosityModel_ == LVISC_INTERACTION) { - m_viscmix = dot_product(m_viscSpecies, m_molefracs); - //now sum over i,j : Gij*Xi*Xj - } else if (viscosityModel_ == LVISC_WILKES) { - err("Wilkes method not implemented for LiquidTransport."); + + // 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); + for ( int i = 0; i < m_nsp; i++ ) + for ( int j = 0; j < i; j++ ) + interaction += m_molefracs[i] * m_molefracs[j] + * ( m_visc_Sij(i,j) + m_visc_Eij(i,j) / m_temp ); + m_viscmix = exp( interaction ); + } return m_viscmix; @@ -540,10 +546,6 @@ namespace Cantera { m_temp = t; m_logt = log(m_temp); m_kbt = Boltzmann * m_temp; - m_sqrt_t = sqrt(m_temp); - m_t14 = sqrt(m_sqrt_t); - m_t32 = m_temp * m_sqrt_t; - m_sqrt_kbt = sqrt(Boltzmann*m_temp); // temperature has changed so temp flags are flipped m_visc_temp_ok = false; @@ -749,9 +751,14 @@ namespace Cantera { int k; for (k = 0; k < m_nsp; k++) { - m_viscSpecies[k] = m_visc_A[k] * exp( m_visc_n[k] * m_logt - - m_visc_Tact[k] / m_temp ); + m_logViscSpecies[k] = m_visc_logA[k] + m_visc_n[k] * m_logt + + m_visc_Tact[k] / m_temp ; + m_viscSpecies[k] = exp( m_logViscSpecies[k] ); } + //for (k = 0; k < m_nsp; k++) { + //m_viscSpecies[k] = m_visc_A[k] * exp( m_visc_n[k] * m_logt + // + m_visc_Tact[k] / m_temp ); + //} m_visc_temp_ok = true; m_visc_mix_ok = false; } diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h index 424338fbe..06f93bb37 100644 --- a/Cantera/src/transport/LiquidTransport.h +++ b/Cantera/src/transport/LiquidTransport.h @@ -31,8 +31,8 @@ using namespace std; namespace Cantera { const int LVISC_CONSTANT = 0; - const int LVISC_WILKES = 1; - const int LVISC_MIXTUREAVG = 2; + const int LVISC_INTERACTION = 1; + const int LVISC_AVG_ENERGIES = 2; const int LDIFF_MIXDIFF_UNCORRECTED = 0; const int LDIFF_MIXDIFF_FLUXCORRECTED = 1; @@ -136,8 +136,13 @@ namespace Cantera { class LiquidTransport : public Transport { public: - //! default constructor + //! Default constructor. /*! + * This requires call to initLiquid(LiquidTransportParams& tr) + * after filling LiquidTransportParams to complete instantiation. + * The filling of LiquidTransportParams is currently carried out + * in the TransportFactory class, but might be moved at some point. + * * @param thermo ThermoPhase object holding species information. * @param ndim Number of spatial dimensions. */ @@ -371,14 +376,32 @@ namespace Cantera { //! Pure species viscosities in Arrhenius temperature-dependent form. vector_fp m_visc_A; + vector_fp m_visc_logA; //logarithm of coefficient vector_fp m_visc_n; vector_fp m_visc_Tact; + //! Molecular interaction energies associated with viscosity + /** + * These multiply the viscosity according to + * \f[ exp( \sum_{i} \sum_{j} X_i X_j E_{i,j} / T \f]. + */ + DenseMatrix m_visc_Eij; + + //! Molecular interaction entropies associated with viscosity + /** + * These multiply the viscosity according to + * \f[ exp( \sum_{i} \sum{j} X_i X_j S_{i,j} \f]. + */ + DenseMatrix m_visc_Sij; + //! Pure species thermal conductivities in Arrhenius temperature-dependent form. vector_fp m_thermCond_A; vector_fp m_thermCond_n; vector_fp m_thermCond_Tact; + //! Species hydrodynamic radius + vector_fp m_hydrodynamic_radius; + //! Polynomial coefficients of the binary diffusion coefficients /*! @@ -469,9 +492,9 @@ namespace Cantera { */ DenseMatrix m_bdiff; - //! Species viscosities + //! Species viscosities and their logarithm /*! - * Viscosity of the species + * Viscosity of the species and its logarithm * Length = number of species * * Depends on the temperature. We have set the pressure dependence @@ -480,6 +503,7 @@ namespace Cantera { * controlling update boolean -> m_visc_temp_ok */ vector_fp m_viscSpecies; + vector_fp m_logViscSpecies; //! Internal value of the species individual thermal conductivities /*! @@ -580,10 +604,6 @@ namespace Cantera { //! Matrix for the stefan maxwell equation. DenseMatrix m_A; - //! Internal storage for the species LJ well depth - vector_fp m_eps; - - //! Current Temperature -> locally storred /*! * This is used to test whether new temperature computations @@ -597,21 +617,6 @@ namespace Cantera { //! Current value of kT doublereal m_kbt; - //! Current Temperature **0.5 - doublereal m_sqrt_t; - - //! Current Temperature **0.25 - doublereal m_t14; - - //! Current Temperature **1.5 - doublereal m_t32; - - //! Current temperature function - /*! - * This is equal to sqrt(Boltzmann * T) - */ - doublereal m_sqrt_kbt; - //! Current value of the pressure doublereal m_press; diff --git a/Cantera/src/transport/LiquidTransportParams.h b/Cantera/src/transport/LiquidTransportParams.h index c3021080f..08142b745 100644 --- a/Cantera/src/transport/LiquidTransportParams.h +++ b/Cantera/src/transport/LiquidTransportParams.h @@ -27,14 +27,37 @@ namespace Cantera { //section for liquid transport properties - //Arrhenius parameters for transport coefficients - // std::vector viscParams; + //Arrhenius parameters for transport coefficients: + + //!Arrhenius pre-exponential parameter for viscosity. vector_fp visc_A; + //!Temperature exponent for viscosity. vector_fp visc_n; + //!Arrhenius activation temperature for viscosity. vector_fp visc_Tact; + + //!Arrhenius pre-exponential parameter for thermal conductivity. vector_fp thermCond_A; + //!Temperature exponent for thermal conductivity. vector_fp thermCond_n; + //!Arrhenius activation temperature for thermal conductivity. vector_fp thermCond_Tact; + + //! Energies of molecular interaction associated with viscosity. + /** + * These multiply the mixture viscosity by + * \f[ \exp( \sum_{i} \sum_{j} X_i X_j ( S_{i,j} + E_{i,j} / T ) ) \f]. + * + * The overall formula for the logarithm of the mixture viscosity is + * + * \f[ \ln \eta_{mix} = \sum_i X_i \ln \eta_i + * + \sum_i \sum_j X_i X_j ( S_{i,j} + E_{i,j} / T ) \f]. + */ + DenseMatrix visc_Eij; + + //! Entropies of molecular interaction associated with viscosity. + DenseMatrix visc_Sij; + //Hydrodynamic radius of transported molecule vector_fp hydroRadius; diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index 19d244fe3..aa7490c2b 100755 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -940,6 +940,8 @@ namespace Cantera { XML_Node& thermCond = trNode.child("thermal_conductivity"); getArrhenius(thermCond, A_thcond, n_thcond, Tact_thcond ); + // Fill datatable with LiquidTransportData objects for error checking + // and then insertion into LiquidTransportData objects below. LiquidTransportData data; data.speciesName = name; @@ -996,6 +998,11 @@ namespace Cantera { trParam.hydroRadius[i] = 1.e-10 * trdat.hydroradius; } + + // Need to identify a method to obtain interaction matrices. + // This will fill LiquidTransportParams members visc_Eij, visc_Sij + trParam.visc_Eij.resize(trParam.nsp_,trParam.nsp_); + cout << "No support for species viscosity interactions in TransportFactory.cpp" << endl; }