From 60afcc5d573df8c0d380f2a4736e7017aeb817b9 Mon Sep 17 00:00:00 2001 From: Christopher Lueth Date: Mon, 1 Mar 2010 16:16:15 +0000 Subject: [PATCH] Implemented stefan maxwell equations with input of conductivity, cation mobility ratio, and self diffusion coefficients. --- Cantera/src/transport/LiquidTransport.cpp | 617 ++++++++++++++-- Cantera/src/transport/LiquidTransport.h | 280 ++++++- Cantera/src/transport/LiquidTransportData.cpp | 100 ++- Cantera/src/transport/LiquidTransportData.h | 78 +- .../src/transport/LiquidTransportParams.cpp | 689 ++++++++++++++---- Cantera/src/transport/LiquidTransportParams.h | 272 +++++-- Cantera/src/transport/TransportBase.h | 70 +- Cantera/src/transport/TransportFactory.cpp | 138 +++- 8 files changed, 1930 insertions(+), 314 deletions(-) diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index ee4756a03..262b12f0f 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -34,23 +34,38 @@ namespace Cantera { LiquidTransport::LiquidTransport(thermo_t* thermo, int ndim) : Transport(thermo, ndim), m_nsp(0), + m_nBinInt(0), m_tmin(-1.0), m_tmax(100000.), m_iStateMF(-1), m_temp(-1.0), m_press(-1.0), + m_lambda(-1.0), m_viscmix(-1.0), + m_ionCondmix(-1.0), + m_mobRatMix(-1.0), + m_selfDiffMix(-1.0), + m_visc_mix_ok(false), m_visc_temp_ok(false), m_visc_conc_ok(false), + m_ionCond_mix_ok(false), + m_ionCond_temp_ok(false), + m_ionCond_conc_ok(false), + m_mobRat_mix_ok(false), + m_mobRat_temp_ok(false), + m_mobRat_conc_ok(false), + m_selfDiff_mix_ok(false), + m_selfDiff_temp_ok(false), + m_selfDiff_conc_ok(false), m_radi_mix_ok(false), m_radi_temp_ok(false), m_radi_conc_ok(false), m_diff_mix_ok(false), m_diff_temp_ok(false), - m_cond_temp_ok(false), - m_cond_mix_ok(false), + m_lambda_temp_ok(false), + m_lambda_mix_ok(false), m_mode(-1000), m_debug(false), m_nDim(1) @@ -61,23 +76,38 @@ namespace Cantera { LiquidTransport::LiquidTransport(const LiquidTransport &right) : Transport(), m_nsp(0), + m_nBinInt(0), m_tmin(-1.0), m_tmax(100000.), m_iStateMF(-1), m_temp(-1.0), m_press(-1.0), + m_lambda(-1.0), m_viscmix(-1.0), + m_ionCondmix(-1.0), + m_mobRatMix(-1.0), + m_selfDiffMix(-1.0), + m_visc_mix_ok(false), m_visc_temp_ok(false), m_visc_conc_ok(false), + m_ionCond_mix_ok(false), + m_ionCond_temp_ok(false), + m_ionCond_conc_ok(false), + m_mobRat_mix_ok(false), + m_mobRat_temp_ok(false), + m_mobRat_conc_ok(false), + m_selfDiff_mix_ok(false), + m_selfDiff_temp_ok(false), + m_selfDiff_conc_ok(false), m_radi_mix_ok(false), m_radi_temp_ok(false), m_radi_conc_ok(false), m_diff_mix_ok(false), m_diff_temp_ok(false), - m_cond_temp_ok(false), - m_cond_mix_ok(false), + m_lambda_temp_ok(false), + m_lambda_mix_ok(false), m_mode(-1000), m_debug(false), m_nDim(1) @@ -95,10 +125,16 @@ namespace Cantera { } Transport::operator=(right); m_nsp = right.m_nsp; + m_nBinInt = right.m_nBinInt; m_tmin = right.m_tmin; m_tmax = right.m_tmax; m_mw = right.m_mw; m_viscTempDep_Ns = right.m_viscTempDep_Ns; + m_ionCondTempDep_Ns = right.m_ionCondTempDep_Ns; + m_mobRatTempDep_Ns = right.m_mobRatTempDep_Ns; + m_mobRatTempDepIndex = right.m_mobRatTempDepIndex; + m_selfDiffTempDep_Ns = right.m_selfDiffTempDep_Ns; + m_selfDiffTempDepIndex = right.m_selfDiffTempDepIndex; m_lambdaTempDep_Ns = right.m_lambdaTempDep_Ns; m_diffTempDep_Ns = right.m_diffTempDep_Ns; m_radiusTempDep_Ns = right.m_radiusTempDep_Ns; @@ -109,9 +145,19 @@ namespace Cantera { m_Grad_mu = right.m_Grad_mu; m_bdiff = right.m_bdiff; m_viscSpecies = right.m_viscSpecies; + m_ionCondSpecies = right.m_ionCondSpecies; + m_mobRatSpecies = right.m_mobRatSpecies; + m_mobRatSpeciesIndex = right.m_mobRatSpeciesIndex; + m_selfDiffSpecies = right.m_selfDiffSpecies; + m_selfDiffSpeciesIndex = right.m_selfDiffSpeciesIndex; m_hydrodynamic_radius = right.m_hydrodynamic_radius; m_lambdaSpecies = right.m_lambdaSpecies; m_viscMixModel = right.m_viscMixModel; + m_ionCondMixModel = right.m_ionCondMixModel; + m_mobRatMixModel = right.m_mobRatMixModel; + m_mobRatMixModelIndex = right.m_mobRatMixModelIndex; + m_selfDiffMixModel = right.m_selfDiffMixModel; + m_selfDiffMixModelIndex = right.m_selfDiffMixModelIndex; m_lambdaMixModel = right.m_lambdaMixModel; m_diffMixModel = right.m_diffMixModel; m_iStateMF = -1; @@ -123,7 +169,6 @@ namespace Cantera { m_actCoeff = right.m_actCoeff; m_Grad_lnAC = right.m_Grad_lnAC; m_chargeSpecies = right.m_chargeSpecies; - m_volume_spec = right.m_volume_spec; m_B = right.m_B; m_A = right.m_A; m_temp = right.m_temp; @@ -132,17 +177,31 @@ namespace Cantera { m_Vdiff = right.m_Vdiff; m_lambda = right.m_lambda; m_viscmix = right.m_viscmix; + m_ionCondmix = right.m_ionCondmix; + m_mobRatMix = right.m_mobRatMix; + m_mobRatMixIndex = right.m_mobRatMixIndex; + m_selfDiffMix = right.m_selfDiffMix; + m_selfDiffMixIndex = right.m_selfDiffMixIndex; m_spwork = right.m_spwork; m_visc_mix_ok = false; m_visc_temp_ok = false; m_visc_conc_ok = false; + m_ionCond_mix_ok = false; + m_ionCond_temp_ok = false; + m_ionCond_conc_ok = false; + m_mobRat_mix_ok = false; + m_mobRat_temp_ok = false; + m_mobRat_conc_ok = false; + m_selfDiff_mix_ok = false; + m_selfDiff_temp_ok = false; + m_selfDiff_conc_ok = false; m_radi_mix_ok = false; m_radi_temp_ok = false; m_radi_conc_ok = false; m_diff_mix_ok = false; m_diff_temp_ok = false; - m_cond_temp_ok = false; - m_cond_mix_ok = false; + m_lambda_temp_ok = false; + m_lambda_mix_ok = false; m_mode = right.m_mode; m_debug = right.m_debug; m_nDim = right.m_nDim; @@ -158,18 +217,33 @@ namespace Cantera { LiquidTransport::~LiquidTransport() { - //These are constructed in TransportFactory::newLTP - for ( int k = 0; k < m_nsp; k++) { - if ( m_viscTempDep_Ns[k] ) delete m_viscTempDep_Ns[k]; - if ( m_lambdaTempDep_Ns[k] ) delete m_lambdaTempDep_Ns[k]; - if ( m_radiusTempDep_Ns[k] ) delete m_radiusTempDep_Ns[k]; - if ( m_diffTempDep_Ns[k] ) delete m_diffTempDep_Ns[k]; - } - //These are constructed in TransportFactory::newLTI - if ( m_viscMixModel ) delete m_viscMixModel; - if ( m_lambdaMixModel ) delete m_lambdaMixModel; - if ( m_diffMixModel ) delete m_diffMixModel; - //if ( m_radiusMixModel ) delete m_radiusMixModel; + //These are constructed in TransportFactory::newLTP + for ( int k = 0; k < m_nsp; k++) { + if ( m_viscTempDep_Ns[k] ) delete m_viscTempDep_Ns[k]; + if ( m_ionCondTempDep_Ns[k] ) delete m_ionCondTempDep_Ns[k]; + for ( int l=0;l < m_nsp; l++ ){ + if ( m_selfDiffTempDep_Ns[l][k] ) delete m_selfDiffTempDep_Ns[l][k]; + } + for ( int l=0;l < m_nBinInt; l++ ){ + if ( m_mobRatTempDep_Ns[l][k] ) delete m_mobRatTempDep_Ns[l][k]; + } + if ( m_lambdaTempDep_Ns[k] ) delete m_lambdaTempDep_Ns[k]; + if ( m_radiusTempDep_Ns[k] ) delete m_radiusTempDep_Ns[k]; + if ( m_diffTempDep_Ns[k] ) delete m_diffTempDep_Ns[k]; + //These are constructed in TransportFactory::newLTI + + if ( m_selfDiffMixModel[k] ) delete m_selfDiffMixModel[k]; + } + + for ( int k = 0; k < m_nBinInt; k++) { + if ( m_mobRatMixModel[k] ) delete m_mobRatMixModel[k]; + } + + if ( m_viscMixModel ) delete m_viscMixModel; + if ( m_ionCondMixModel ) delete m_ionCondMixModel; + if ( m_lambdaMixModel ) delete m_lambdaMixModel; + if ( m_diffMixModel ) delete m_diffMixModel; + //if ( m_radiusMixModel ) delete m_radiusMixModel; } @@ -191,6 +265,7 @@ namespace Cantera { m_thermo = tr.thermo; m_velocityBasis = tr.velocityBasis_; m_nsp = m_thermo->nSpecies(); + m_nBinInt = m_nsp*(m_nsp-1)/2; m_tmin = m_thermo->minTemp(); m_tmax = m_thermo->maxTemp(); @@ -200,38 +275,77 @@ namespace Cantera { m_thermo->molecularWeights().end(), m_mw.begin()); /* - * Get the input Viscosities + * Get the input Viscosities, and stuff */ m_viscSpecies.resize(m_nsp, 0.0); m_viscTempDep_Ns.resize(m_nsp, 0); + m_ionCondSpecies.resize(m_nsp, 0.0); + m_ionCondTempDep_Ns.resize(m_nsp, 0); + m_mobRatTempDepIndex.resize(m_nBinInt); + m_mobRatTempDep_Ns.resize(m_nBinInt); + m_mobRatMixModel.resize(m_nBinInt); + m_mobRatMixModelIndex.resize(m_nBinInt); + m_mobRatSpeciesIndex.resize(m_nBinInt); + m_mobRatSpecies.resize(m_nBinInt, m_nsp, 0.0); + m_mobRatMix.resize(m_nBinInt,0.0); + m_mobRatMixIndex.resize(m_nBinInt); + m_selfDiffTempDepIndex.resize(m_nsp); + m_selfDiffTempDep_Ns.resize(m_nsp); + m_selfDiffMixModel.resize(m_nsp); + m_selfDiffMixModelIndex.resize(m_nsp); + m_selfDiffSpeciesIndex.resize(m_nsp); + m_selfDiffSpecies.resize(m_nsp, m_nsp, 0.0); + m_selfDiffMix.resize(m_nsp,0.0); + m_selfDiffMixIndex.resize(m_nsp); + for (k=0; k < m_nsp; k++){ + m_selfDiffTempDep_Ns[k].resize(m_nsp, 0); + } + for (k=0; k < m_nBinInt; k++){ + m_mobRatTempDep_Ns[k].resize(m_nsp, 0); + } + m_lambdaSpecies.resize(m_nsp, 0.0); + m_lambdaTempDep_Ns.resize(m_nsp, 0); + m_hydrodynamic_radius.resize(m_nsp, 0.0); + m_radiusTempDep_Ns.resize(m_nsp, 0); + + //first populate mixing rules and indices + for (k = 0; k < m_nsp; k++) { + m_selfDiffMixModel[k] = tr.selfDiffusion[k]; + m_selfDiffMixModelIndex[k] = tr.selfDiffIndex[k]; + } + for (k = 0; k < m_nBinInt; k++) { + m_mobRatMixModel[k] = tr.mobilityRatio[k]; + m_mobRatMixModelIndex[k] = tr.mobRatIndex[k]; + } + //for each species, assign viscosity model and coefficients for (k = 0; k < m_nsp; k++) { Cantera::LiquidTransportData <d = tr.LTData[k]; m_viscTempDep_Ns[k] = ltd.viscosity; - } - - /* - * Get the input Thermal Conductivities - */ - m_lambdaSpecies.resize(m_nsp, 0.0); - m_lambdaTempDep_Ns.resize(m_nsp, 0); - //for each species, assign thermal conductivity model - for (k = 0; k < m_nsp; k++) { - Cantera::LiquidTransportData <d = tr.LTData[k]; + m_ionCondTempDep_Ns[k] = ltd.ionConductivity; + for (int j = 0; j < m_nBinInt; j++){ + for (int l=0; l < m_nBinInt; l++){ + if (m_mobRatMixModelIndex[j] == ltd.mobRatIndex[l]) { + m_mobRatTempDep_Ns[j][k] = ltd.mobilityRatio[l]; + m_mobRatTempDepIndex[j] = ltd.mobRatIndex[l]; + break; + } + } + } + for (int j = 0; j < (int) m_selfDiffMixModelIndex.size(); j++){ + for (int l=0; l < (int) m_selfDiffMixModelIndex.size(); l++){ + if (m_selfDiffMixModelIndex[j] == ltd.selfDiffIndex[l]) { + m_selfDiffTempDep_Ns[j][k] = ltd.selfDiffusion[l]; + m_selfDiffTempDepIndex[j] = ltd.selfDiffIndex[l]; + break; + } + } + } m_lambdaTempDep_Ns[k] = ltd.thermalCond; - } - - /* - * Get the input Hydrodynamic Radii - */ - m_hydrodynamic_radius.resize(m_nsp, 0.0); - m_radiusTempDep_Ns.resize(m_nsp, 0); - //for each species, assign model for hydrodynamic radius - for (k = 0; k < m_nsp; k++) { - Cantera::LiquidTransportData <d = tr.LTData[k]; m_radiusTempDep_Ns[k] = ltd.hydroRadius; } + /* * Get the input Species Diffusivities * Note that species diffusivities are not what is needed. @@ -263,18 +377,28 @@ namespace Cantera { * species diffusivity and hydrodynamics radius (perhaps not needed in the * present class). */ + + m_viscMixModel = tr.viscosity; tr.viscosity = 0; + + m_ionCondMixModel = tr.ionConductivity; + tr.ionConductivity = 0; + //m_mobRatMixModel = tr.mobilityRatio; + m_lambdaMixModel = tr.thermalCond; tr.thermalCond = 0; m_diffMixModel = tr.speciesDiffusivity; tr.speciesDiffusivity = 0; - + + m_bdiff.resize(m_nsp,m_nsp, 0.0); + //Don't really need to update this here. //It is updated in updateDiff_T() - m_diffMixModel->getMatrixTransProp(m_bdiff); + m_diffMixModel->getMatrixTransProp( m_bdiff ); + m_mode = tr.mode_; @@ -306,10 +430,19 @@ namespace Cantera { m_visc_mix_ok = false; m_visc_temp_ok = false; m_visc_conc_ok = false; + m_ionCond_mix_ok = false; + m_ionCond_temp_ok = false; + m_ionCond_conc_ok = false; + m_mobRat_mix_ok = false; + m_mobRat_temp_ok = false; + m_mobRat_conc_ok = false; + m_selfDiff_mix_ok = false; + m_selfDiff_temp_ok = false; + m_selfDiff_conc_ok = false; m_radi_temp_ok = false; m_radi_conc_ok = false; - m_cond_temp_ok = false; - m_cond_mix_ok = false; + m_lambda_temp_ok = false; + m_lambda_mix_ok = false; m_diff_temp_ok = false; m_diff_mix_ok = false; @@ -357,6 +490,215 @@ namespace Cantera { copy(m_viscSpecies.begin(), m_viscSpecies.end(), visc); } + /****************** ionConductivity ******************************/ + + // Returns the ionic conductivity of the solution + /* + * The ionConductivity calculation is handled by subclasses of + * LiquidTranInteraction as specified in the input file. + * These in turn employ subclasses of LTPspecies to + * determine the individual species ionic conductivities. + */ + + doublereal LiquidTransport:: ionConductivity() { + + update_T(); + update_C(); + + if (m_ionCond_mix_ok) return m_ionCondmix; + + ////// LiquidTranInteraction method + m_ionCondmix = m_ionCondMixModel->getMixTransProp( m_ionCondTempDep_Ns ); + + return m_ionCondmix; + + /* + // update m_ionCondSpecies[] if necessary + if (!m_ionCond_temp_ok) { + updateIonConductivity_T(); + } + + if (!m_ionCond_conc_ok) { + updateIonConductivity_C(); + } + */ + } + + // Returns the pure species ionic conductivities for all species + /* + * The pure species ionic conductivities are evaluated using the + * appropriate subclasses of LTPspecies as specified in the + * input file. + * + * @param ionCond array of length "number of species" + * to hold returned ionic conductivities. + */ + void LiquidTransport::getSpeciesIonConductivity(doublereal* ionCond) { + update_T(); + if (!m_ionCond_temp_ok) { + updateIonConductivity_T(); + } + copy(m_ionCondSpecies.begin(), m_ionCondSpecies.end(), ionCond); + } + + /****************** mobilityRatio ******************************/ + + // Returns the mobility ratios of the solution + /* + * The mobility ratio calculation is handled by subclasses of + * LiquidTranInteraction as specified in the input file. + * These in turn employ subclasses of LTPspecies to + * determine the individual species mobility ratios. + */ + void LiquidTransport:: mobilityRatio(vector_fp& mobRat, std::vector& mobRatIndex) { + + update_T(); + update_C(); + + ////// LiquidTranInteraction method + if (!m_mobRat_mix_ok){ + for (int k = 0; k < m_nBinInt; k++){ + if ( m_mobRatMixModelIndex[k] != m_mobRatTempDepIndex[k] ) + throw CanteraError("LiquidTransport::mobilityRation","Mobility Ratio Indices Don't Match: Mixture vs. Species"); + m_mobRatMix[k] = m_mobRatMixModel[k]->getMixTransProp( m_mobRatTempDep_Ns[k] ); + m_mobRatMixIndex[k] = m_mobRatMixModelIndex[k]; + } + } + for (int k = 0; k < m_nBinInt; k++){ + mobRat[k] = m_mobRatMix[k]; + mobRatIndex[k]= m_mobRatMixIndex[k]; + } + } + void LiquidTransport:: mobilityRatio(doublereal* mobRat, std::vector& mobRatIndex) { + + update_T(); + update_C(); + + ////// LiquidTranInteraction method + if (!m_mobRat_mix_ok){ + for (int k = 0; k < m_nBinInt; k++){ + if ( m_mobRatMixModelIndex[k] != m_mobRatTempDepIndex[k] ) + throw CanteraError("LiquidTransport::mobilityRatio","Mobility Ratio Indices Don't Match: Mixture vs. Species"); + m_mobRatMix[k] = m_mobRatMixModel[k]->getMixTransProp( m_mobRatTempDep_Ns[k] ); + m_mobRatMixIndex[k] = m_mobRatMixModelIndex[k]; + } + } + for (int k = 0; k < m_nBinInt; k++){ + mobRat[k] = m_mobRatMix[k]; + mobRatIndex[k]= m_mobRatMixIndex[k]; + } + } + + // Returns the pure species mobility ratios for all species + /* + * The pure species mobility ratios are evaluated using the + * appropriate subclasses of LTPspecies as specified in the + * input file. + * + * @param mobRat array of length "number of species" + * to hold returned mobility ratio. + */ + void LiquidTransport::getSpeciesMobilityRatio(DenseMatrix& mobRat, std::vector& mobRatIndex) { + update_T(); + if (!m_mobRat_temp_ok) { + updateMobilityRatio_T(); + } + mobRat = m_mobRatSpecies; + for (int k = 0; k < m_nBinInt; k++) + mobRatIndex[k] = m_mobRatSpeciesIndex[k]; + } + void LiquidTransport::getSpeciesMobilityRatio(doublereal** mobRat, std::vector& mobRatIndex) { + update_T(); + if (!m_mobRat_temp_ok) { + updateMobilityRatio_T(); + } + for (int k=0; k& selfDiffIndex) { + + update_T(); + update_C(); + + ////// LiquidTranInteraction method + if (!m_selfDiff_mix_ok){ + for (int k = 0; k < m_nsp; k++){ + if ( m_selfDiffMixModelIndex[k] != m_selfDiffTempDepIndex[k] ) + throw CanteraError("LiquidTransport::selfDiffusion","Self Diffusion Indices Don't Match: Mixture vs. Species"); + m_selfDiffMix[k] = m_selfDiffMixModel[k]->getMixTransProp( m_selfDiffTempDep_Ns[k] ); + m_selfDiffMixIndex[k] = m_selfDiffMixModelIndex[k]; + } + } + for (int k = 0; k < m_nsp; k++){ + selfDiff[k] = m_selfDiffMix[k]; + selfDiffIndex[k]= m_selfDiffMixIndex[k]; + } + } + + void LiquidTransport:: selfDiffusion(double* selfDiff, std::vector& selfDiffIndex) { + + update_T(); + update_C(); + + ////// LiquidTranInteraction method + if (!m_selfDiff_mix_ok){ + for (int k = 0; k < m_nsp; k++){ + if ( m_selfDiffMixModelIndex[k] != m_selfDiffTempDepIndex[k] ) + throw CanteraError("LiquidTransport::selfDiffusion","Self Diffusion Indices Don't Match: Mixture vs. Species"); + m_selfDiffMix[k] = m_selfDiffMixModel[k]->getMixTransProp( m_selfDiffTempDep_Ns[k] ); + m_selfDiffMixIndex[k] = m_selfDiffMixModelIndex[k]; + } + } + for (int k = 0; k < m_nsp; k++){ + selfDiff[k] = m_selfDiffMix[k]; + selfDiffIndex[k]= m_selfDiffMixIndex[k]; + } + } + + // Returns the pure species self diffusion for all species + /* + * The pure species self diffusion coeffs are evaluated using the + * appropriate subclasses of LTPspecies as specified in the + * input file. + * + * @param selfDiff array of size "number of species"^2 + * to hold returned self diffusion. + */ + void LiquidTransport::getSpeciesSelfDiffusion(DenseMatrix& selfDiff, std::vector& selfDiffIndex) { + update_T(); + if (!m_selfDiff_temp_ok) { + updateSelfDiffusion_T(); + } + selfDiff = m_selfDiffSpecies; + for (int k = 0; k < m_nsp; k++) + selfDiffIndex[k] = m_selfDiffSpeciesIndex[k]; + } + void LiquidTransport::getSpeciesSelfDiffusion(doublereal** selfDiff, std::vector& selfDiffIndex) { + update_T(); + if (!m_selfDiff_temp_ok) { + updateSelfDiffusion_T(); + } + for (int k=0; kgetMixTransProp( m_lambdaTempDep_Ns ); m_cond_mix_ok = true; } @@ -846,11 +1188,6 @@ namespace Cantera { void LiquidTransport::getSpeciesVdiffExt(int ldf, doublereal* Vdiff) { int n, k; - update_T(); - update_C(); - - update_Grad_lnAC(); - stefan_maxwell_solve(); for (n = 0; n < m_nDim; n++) { @@ -878,11 +1215,6 @@ namespace Cantera { void LiquidTransport::getSpeciesFluxesExt(int ldf, doublereal* fluxes) { int n, k; - update_T(); - update_C(); - - update_Grad_lnAC(); - stefan_maxwell_solve(); for (n = 0; n < m_nDim; n++) { @@ -926,11 +1258,6 @@ namespace Cantera { */ void LiquidTransport::getMixDiffCoeffs(doublereal* const d) { - update_T(); - update_C(); - - update_Grad_lnAC(); - stefan_maxwell_solve(); for ( int n = 0; n < m_nDim; n++) { @@ -975,19 +1302,28 @@ namespace Cantera { // temperature has changed so temp flags are flipped m_visc_temp_ok = false; + m_ionCond_temp_ok = false; + m_mobRat_temp_ok = false; + m_selfDiff_temp_ok = false; m_radi_temp_ok = false; m_diff_temp_ok = false; + m_lambda_temp_ok = false; // temperature has changed, so polynomial temperature // interpolations will need to be reevaluated. // This means that many concentration m_visc_conc_ok = false; - m_cond_temp_ok = false; + m_ionCond_conc_ok = false; + m_mobRat_conc_ok = false; + m_selfDiff_conc_ok = false; // Mixture stuff needs to be evaluated m_visc_mix_ok = false; + m_ionCond_mix_ok = false; + m_mobRat_mix_ok = false; + m_selfDiff_mix_ok = false; m_diff_mix_ok = false; - // m_cond_mix_ok = false; (don't need it because a lower lvl flag is set + m_lambda_mix_ok = false; //(don't need it because a lower lvl flag is set return true; } @@ -1040,11 +1376,17 @@ namespace Cantera { // be recomputed before use, and update the local mole // fractions. m_visc_conc_ok = false; + m_ionCond_conc_ok = false; + m_mobRat_conc_ok = false; + m_selfDiff_conc_ok = false; // Mixture stuff needs to be evaluated m_visc_mix_ok = false; + m_ionCond_mix_ok = false; + m_mobRat_mix_ok = false; + m_selfDiff_mix_ok = false; m_diff_mix_ok = false; - m_cond_mix_ok = false; + m_lambda_mix_ok = false; return true; } @@ -1067,8 +1409,8 @@ namespace Cantera { for (k = 0; k < m_nsp; k++) { m_lambdaSpecies[k] = m_lambdaTempDep_Ns[k]->getSpeciesTransProp() ; } - m_cond_temp_ok = true; - m_cond_mix_ok = false; + m_lambda_temp_ok = true; + m_lambda_mix_ok = false; } @@ -1112,6 +1454,79 @@ namespace Cantera { } + //! Update the pure-species ionic conductivities functional dependence on concentration. + void LiquidTransport::updateIonConductivity_C() { + m_ionCond_conc_ok = true; + } + + + /** + * Updates the array of pure species ionic conductivities internally + * using calls to the appropriate LTPspecies subclass. + * The flag m_ionCond_ok is set to true. + */ + void LiquidTransport::updateIonConductivity_T() { + int k; + + for (k = 0; k < m_nsp; k++) { + m_ionCondSpecies[k] = m_ionCondTempDep_Ns[k]->getSpeciesTransProp() ; + } + m_ionCond_temp_ok = true; + m_ionCond_mix_ok = false; + } + + + //! Update the pure-species mobility ratios functional dependence on concentration. + void LiquidTransport::updateMobilityRatio_C() { + m_mobRat_conc_ok = true; + } + + + /** + * Updates the array of pure species mobility ratios internally + * using calls to the appropriate LTPspecies subclass. + * The flag m_mobRat_ok is set to true. + */ + void LiquidTransport::updateMobilityRatio_T() { + int k; + int j; + + for (k = 0; k < m_nBinInt; k++) { + for (j = 0; j < m_nsp; j++) { + m_mobRatSpecies(k,j) = m_mobRatTempDep_Ns[k][j]->getSpeciesTransProp() ; + } + m_mobRatSpeciesIndex[k] = m_mobRatTempDepIndex[k]; + } + m_mobRat_temp_ok = true; + m_mobRat_mix_ok = false; + } + + + //! Update the pure-species self diffusion functional dependence on concentration. + void LiquidTransport::updateSelfDiffusion_C() { + m_selfDiff_conc_ok = true; + } + + + /** + * Updates the array of pure species self diffusion internally + * using calls to the appropriate LTPspecies subclass. + * The flag m_selfDiff_ok is set to true. + */ + void LiquidTransport::updateSelfDiffusion_T() { + int k; + int j; + + for (k = 0; k < m_nBinInt; k++) { + for (j = 0; j < m_nsp; j++) { + m_selfDiffSpecies(k,j) = m_selfDiffTempDep_Ns[k][j]->getSpeciesTransProp() ; + } + m_selfDiffSpeciesIndex[k] = m_selfDiffTempDepIndex[k]; + } + m_selfDiff_temp_ok = true; + m_selfDiff_mix_ok = false; + } + //! Update the pure-species viscosities functional dependence on concentration. void LiquidTransport::updateHydrodynamicRadius_C() { m_radi_conc_ok = true; @@ -1164,14 +1579,23 @@ namespace Cantera { int k; - vector_fp grad_lnAC(m_nsp); - m_thermo->getdlnActCoeffdlnC( DATA_PTR(grad_lnAC) ); + doublereal grad_T; + vector_fp grad_lnAC(m_nsp), grad_X(m_nsp); + // IonsFromNeutralVPSSTP * tempIons = dynamic_cast m_thermo; + //MargulesVPSSTP * tempMarg = dynamic_cast (tempIons->neutralMoleculePhase_); + - for (k = 0; k < m_nsp; k++) { - m_Grad_lnAC[k] = grad_lnAC[k]; + //m_thermo->getdlnActCoeffdlnX( DATA_PTR(grad_lnAC) ); + for (k = 0; k < m_nDim; k++ ) { + grad_T = m_Grad_T[k]; + grad_X.assign(m_Grad_X.begin()+m_nsp*k,m_Grad_X.begin()+m_nsp*(k+1)); + m_thermo->getdlnActCoeff( grad_T, DATA_PTR(grad_X), DATA_PTR(grad_lnAC) ); + for ( int i = 0; i < m_nsp; i++ ) + grad_lnAC[i] += grad_X[i]/m_molefracs[i]; + copy(grad_lnAC.begin(),grad_lnAC.end(),m_Grad_lnAC.begin()+m_nsp*k); // std::cout << k << " m_Grad_lnAC = " << m_Grad_lnAC[k] << std::endl; } - + return; } @@ -1220,10 +1644,13 @@ namespace Cantera { //! grab a local copy of the molecular weights const vector_fp& M = m_thermo->molecularWeights(); + //! grad a local copy of the ion molar volume (inverse total ion concentration) + const doublereal vol = m_thermo->molarVolume(); /* - * Update the concentrations and diffusion coefficients in the mixture. + * Update the temperature, concentrations and diffusion coefficients in the mixture. */ + update_T(); update_C(); if ( !m_diff_temp_ok ) updateDiff_T(); @@ -1260,12 +1687,11 @@ namespace Cantera { * */ for (i = 0; i < m_nsp; i++) { - double xi_denom = m_molefracs_tran[i]; for (a = 0; a < m_nDim; a++) { m_Grad_mu[a*m_nsp + i] = m_chargeSpecies[i] * Faraday * m_Grad_V[a] //+ (m_volume_spec[i] - M[i]/dens_) * m_Grad_P[a] - + GasConstant * T * m_Grad_X[a*m_nsp+i] * ( 1.0 + m_Grad_lnAC[i] ) / xi_denom; + + GasConstant * T * m_Grad_lnAC[a*m_nsp+i]; } } @@ -1286,8 +1712,11 @@ 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. */ + + double condSum1, condSum2; switch (m_nDim) { case 1: /* 1-D approximation */ + m_B(0,0) = 0.0; //equation for the reference velocity for (j = 0; j < m_nsp; j++) { @@ -1298,9 +1727,10 @@ namespace Cantera { else if ( ( m_velocityBasis >= 0 ) && ( m_velocityBasis < m_nsp ) ) // use species number m_velocityBasis as reference velocity - if ( m_velocityBasis == j ) m_A(0,j) = 1.0; - else - throw CanteraError("LiquidTransport::stefan_maxwell_solve", + if ( m_velocityBasis == j ) m_A(0,j) = 1.0; + else m_A(0,j) = 0.0; + else + throw CanteraError("LiquidTransport::stefan_maxwell_solve", "Unknown reference velocity provided."); } for (i = 1; i < m_nsp; i++){ @@ -1320,7 +1750,30 @@ namespace Cantera { //! invert and solve the system Ax = b. Answer is in m_B solve(m_A, m_B); - + + /* + condSum2 = m_chargeSpecies[1]*m_chargeSpecies[1]*m_molefracs_tran[1]/m_bdiff(2,3) + + m_chargeSpecies[2]*m_chargeSpecies[2]*m_molefracs_tran[2]/m_bdiff(1,3) + + m_chargeSpecies[3]*m_chargeSpecies[3]*m_molefracs_tran[3]/m_bdiff(1,2); + condSum1 = m_molefracs_tran[1]/m_bdiff(1,2)/m_bdiff(1,3) + + m_molefracs_tran[2]/m_bdiff(2,3)/m_bdiff(1,2) + + m_molefracs_tran[3]/m_bdiff(1,3)/m_bdiff(2,3); + condSum2 = condSum2/condSum1*Faraday*Faraday/GasConstant/T/vol; + */ + + condSum1 = 0; + for (i = 0; i < m_nsp; i++){ + condSum1 -= Faraday*m_chargeSpecies[i]*m_B(i,0)*m_molefracs_tran[i]/vol; + } + + /* + Check Mobility Ratio of Cations + cout << "mobility ratio = " << m_chargeSpecies[1]*(m_B(1,0)-m_B(2,0))/m_chargeSpecies[0]/(m_B(0,0)-m_B(2,0)) << endl; + */ + + // cout << condSum1 << " = " << condSum2 << endl; + + break; case 2: /* 2-D approximation */ m_B(0,0) = 0.0; @@ -1335,8 +1788,9 @@ namespace Cantera { && ( m_velocityBasis < m_nsp ) ) // use species number m_velocityBasis as reference velocity if ( m_velocityBasis == j ) m_A(0,j) = 1.0; - else - throw CanteraError("LiquidTransport::stefan_maxwell_solve", + else m_A(0,j) = 0.0; + else + throw CanteraError("LiquidTransport::stefan_maxwell_solve", "Unknown reference velocity provided."); } for (i = 1; i < m_nsp; i++){ @@ -1375,8 +1829,9 @@ namespace Cantera { && ( m_velocityBasis < m_nsp ) ) // use species number m_velocityBasis as reference velocity if ( m_velocityBasis == j ) m_A(0,j) = 1.0; - else - throw CanteraError("LiquidTransport::stefan_maxwell_solve", + else m_A(0,j) = 0.0; + else + throw CanteraError("LiquidTransport::stefan_maxwell_solve", "Unknown reference velocity provided."); } for (i = 1; i < m_nsp; i++){ diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h index 063c37d0b..133624ca5 100644 --- a/Cantera/src/transport/LiquidTransport.h +++ b/Cantera/src/transport/LiquidTransport.h @@ -68,7 +68,7 @@ namespace Cantera { * \f] * for mass fraction \f$ Y_i \f$. For mole-averaged velocities * \f[ - * \sum_{i} X_i \vec{V_i} = 0 + * \sum_{i} X_i \vec{V_i} = 0 * \f] * for mole fraction \f$ X_i \f$. * @@ -165,7 +165,7 @@ namespace Cantera { * determine the individual species viscosities. */ virtual doublereal viscosity(); - + //! Returns the pure species viscosities for all species /*! * The pure species viscosities are evaluated using the @@ -177,6 +177,70 @@ namespace Cantera { */ virtual void getSpeciesViscosities(doublereal* const visc); + //! Returns the ionic conductivity of the solution + /*! + * The ionic conductivity calculation is handled by subclasses of + * LiquidTranInteraction as specified in the input file. + * These in turn employ subclasses of LTPspecies to + * determine the individual species ionic conductivities. + */ + virtual doublereal ionConductivity(); + + //! Returns the pure species ionic conductivities for all species + /*! + * The pure species ionic conductivities are evaluated using the + * appropriate subclasses of LTPspecies as specified in the + * input file. + * + * @param visc array of length "number of species" + * to hold returned ionic conductivities. + */ + virtual void getSpeciesIonConductivity(doublereal* const ionCond); + + //! Returns the mobility ratio of the solution + /*! + * The mobility ratio calculation is handled by subclasses of + * LiquidTranInteraction as specified in the input file. + * These in turn employ subclasses of LTPspecies to + * determine the individual species mobility ratios. + */ + virtual void mobilityRatio(vector_fp& mobRat, std::vector& mobRatIndex); + virtual void mobilityRatio(double* mobRat, std::vector& mobRatIndex); + + //! Returns the pure species mobility ratios for all species + /*! + * The pure species mobility ratios are evaluated using the + * appropriate subclasses of LTPspecies as specified in the + * input file. + * + * @param mobRat array of length "number of species" + * to hold returned mobility ratios. + */ + virtual void getSpeciesMobilityRatio(DenseMatrix& mobRat, std::vector& mobRatIndex); + virtual void getSpeciesMobilityRatio(double** mobRat, std::vector& mobRatIndex); + + //! Returns the self diffusion coefficients in the solution + /*! + * The self diffusion calculation is handled by subclasses of + * LiquidTranInteraction as specified in the input file. + * These in turn employ subclasses of LTPspecies to + * determine the individual species self diffusion coeffs. + */ + virtual void selfDiffusion(vector_fp& selfDiff, std::vector& selfDiffIndex); + virtual void selfDiffusion(double* selfDiff, std::vector& selfDiffIndex); + + //! Returns the pure species self diffusion in solution of each species + /*! + * The pure species molar volumes are evaluated using the + * appropriate subclasses of LTPspecies as specified in the + * input file. + * + * @param selfDiff array of length "number of species" + * to hold returned self diffusion coeffs. + */ + virtual void getSpeciesSelfDiffusion(DenseMatrix& selfDiff, std::vector& selfDiffIndex); + virtual void getSpeciesSelfDiffusion(double** selfDiff, std::vector& selfDiffIndex); + //! Returns the hydrodynamic radius for all species /*! * The species hydrodynamic radii are evaluated using the @@ -711,6 +775,25 @@ namespace Cantera { */ void updateViscosity_T(); + //! Update the temperature-dependent ionic conductivity terms + //! for each species internally + /*! + * The flag m_ionCond_temp_ok is set to true. + */ + void updateIonConductivity_T(); + + //! Updates the array of pure species mobility ratios internally. + /*! + * The flag m_mobRat_ok is set to true. + */ + void updateMobilityRatio_T(); + + //! Updates the array of pure species self diffusion coeffs internally. + /*! + * The flag m_selfDiff_ok is set to true. + */ + void updateSelfDiffusion_T(); + //! Update the temperature-dependent hydrodynamic radius terms //! for each species internally /*! @@ -731,7 +814,37 @@ namespace Cantera { * @internal */ void updateViscosities_C(); + + //! Update the concentration parts of the ionic conductivity + /*! + * Internal routine is run whenever the update_boolean + * m_ionCond_conc_ok is false. Currently there is no concentration + * dependence for the pure species ionic conductivity. + * + * @internal + */ + void updateIonConductivity_C(); + //! Update the concentration parts of the mobility ratio + /*! + * Internal routine is run whenever the update_boolean + * m_mobRat_conc_ok is false. Currently there is no concentration + * dependence for the pure species mobility ratio. + * + * @internal + */ + void updateMobilityRatio_C(); + + //! Update the concentration parts of the self diffusion + /*! + * Internal routine is run whenever the update_boolean + * m_selfDiff_conc_ok is false. Currently there is no concentration + * dependence for the pure species self diffusion. + * + * @internal + */ + void updateSelfDiffusion_C(); + //! Update the concentration dependence of the hydrodynamic radius /*! * Internal routine is run whenever the update_boolean @@ -751,6 +864,7 @@ namespace Cantera { //! Number of species in the mixture int m_nsp; + int m_nBinInt; //! Minimum temperature applicable to the transport property eval doublereal m_tmin; @@ -782,6 +896,65 @@ namespace Cantera { */ LiquidTranInteraction *m_viscMixModel; + //! Ionic conductivity for each species expressed as an appropriate subclass + //! of LTPspecies + /*! + * These subclasses of LTPspecies evaluate the species-specific + * transport properties according to the parameters parsed in + * TransportFactory::getLiquidSpeciesTransportData(). + */ + std::vector m_ionCondTempDep_Ns; + + //! Ionic Conductivity of the mixture expressed as a subclass of + //! LiquidTranInteraction + /*! + * These subclasses of LiquidTranInteraction evaluate the + * mixture transport properties according to the parameters parsed in + * TransportFactory::getLiquidInteractionsTransportData(). + */ + LiquidTranInteraction *m_ionCondMixModel; + + //! Mobility ratio for each species expressed as an appropriate subclass + //! of LTPspecies + /*! + * These subclasses of LTPspecies evaluate the species-specific + * transport properties according to the parameters parsed in + * TransportFactory::getLiquidSpeciesTransportData(). + */ + typedef std::vector LTPvector; + std::vector m_mobRatTempDep_Ns; + std::vector m_mobRatTempDepIndex; + + //! Mobility ratio of the mixture expressed as a subclass of + //! LiquidTranInteraction + /*! + * These subclasses of LiquidTranInteraction evaluate the + * mixture transport properties according to the parameters parsed in + * TransportFactory::getLiquidInteractionsTransportData(). + */ + std::vector m_mobRatMixModel; + std::vector m_mobRatMixModelIndex; + + //! Self Diffusion for each species expressed as an appropriate subclass + //! of LTPspecies + /*! + * These subclasses of LTPspecies evaluate the species-specific + * transport properties according to the parameters parsed in + * TransportFactory::getLiquidSpeciesTransportData(). + */ + std::vector m_selfDiffTempDep_Ns; + std::vector m_selfDiffTempDepIndex; + + //! Self Diffusion of the mixture expressed as a subclass of + //! LiquidTranInteraction + /*! + * These subclasses of LiquidTranInteraction evaluate the + * mixture transport properties according to the parameters parsed in + * TransportFactory::getLiquidInteractionsTransportData(). + */ + std::vector m_selfDiffMixModel; + std::vector m_selfDiffMixModelIndex; + //! Thermal conductivity for each species expressed as an //! appropriate subclass of LTPspecies /*! @@ -872,19 +1045,24 @@ 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_lnAC[k] = \partial \left[ \ln ( \gamma_i ) \right] - * / \partial \left[ \ln ( \X_i ) \right] + * m_Grad_lnAC[k] = \nabla \ln ( \gamma_i ) + \nabla \ln ( \X_i ) * * Note that where "mole fraction" 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. + * m_nsp is the number of species in the fluid + * + * k is the species index + * n is the dimensional index (x, y, or z). It has a length + * equal to m_nDimm_ + * + * m_Grad_X[n*m_nsp + k] * */ vector_fp m_Grad_lnAC; @@ -964,6 +1142,41 @@ namespace Cantera { */ vector_fp m_viscSpecies; + //! Internal value of the species ionic conductivities + /*! + * Ionic conductivity of the species evaluated using subclass of LTPspecies + * held in m_ionCondTempDep_Ns. + * + * Length = number of species + * + * controlling update boolean -> m_ionCond_temp_ok + */ + vector_fp m_ionCondSpecies; + + //! Internal value of the species mobility ratios + /*! + * Mobility ratio of the species evaluated using subclass of LTPspecies + * held in m_mobRatTempDep_Ns. + * + * Length = number of species + * + * controlling update boolean -> m_mobRat_temp_ok + */ + DenseMatrix m_mobRatSpecies; + std::vector m_mobRatSpeciesIndex; + + //! Internal value of the species self diffusion coefficients + /*! + * Self diffusion of the species evaluated using subclass of LTPspecies + * held in m_selfDiffTempDep_Ns. + * + * Length = number of species + * + * controlling update boolean -> m_selfDiff_temp_ok + */ + DenseMatrix m_selfDiffSpecies; + std::vector m_selfDiffSpeciesIndex; + //! Internal value of the species individual thermal conductivities /*! * Thermal conductivities of the species evaluated using subclass @@ -1093,6 +1306,19 @@ namespace Cantera { //! Saved value of the mixture viscosity doublereal m_viscmix; + //! Saved value of the mixture ionic conductivity + doublereal m_ionCondmix; + + //! Saved values of the mixture mobility ratios + vector_fp m_mobRatMix; + //! Saved species index of the mixture correlated to mobility ratios + std::vector m_mobRatMixIndex; + + //! Saved values of the mixture self diffusion coefficients + vector_fp m_selfDiffMix; + //! Saved species index of the mixture correlated to self diffusion coefficients + std::vector m_selfDiffMixIndex; + //! work space /*! * Length is equal to m_nsp @@ -1115,6 +1341,46 @@ namespace Cantera { //! are current wrt the concentration bool m_visc_conc_ok; + //! Boolean indicating that the top-level mixture ionic conductivity is current + /*! + * This is turned false for every change in T, P, or C. + */ + bool m_ionCond_mix_ok; + + //! Boolean indicating that weight factors wrt ionic conductivty is current + bool m_ionCond_temp_ok; + + //! Flag to indicate that the pure species ionic conductivities + //! are current wrt the concentration + bool m_ionCond_conc_ok; + bool m_cond_mix_ok; + + //! Boolean indicating that the top-level mixture mobility ratio is current + /*! + * This is turned false for every change in T, P, or C. + */ + bool m_mobRat_mix_ok; + + //! Boolean indicating that weight factors wrt mobility ratio is current + bool m_mobRat_temp_ok; + + //! Flag to indicate that the pure species mobility ratios + //! are current wrt the concentration + bool m_mobRat_conc_ok; + + //! Boolean indicating that the top-level mixture self diffusion is current + /*! + * This is turned false for every change in T, P, or C. + */ + bool m_selfDiff_mix_ok; + + //! Boolean indicating that weight factors wrt self diffusion is current + bool m_selfDiff_temp_ok; + + //! Flag to indicate that the pure species self diffusion + //! are current wrt the concentration + bool m_selfDiff_conc_ok; + //! Boolean indicating that mixture diffusion coeffs are current bool m_radi_mix_ok; @@ -1134,10 +1400,10 @@ namespace Cantera { //! Flag to indicate that the pure species conductivities //! are current wrt the temperature - bool m_cond_temp_ok; + bool m_lambda_temp_ok; //! Boolean indicating that mixture conductivity is current - bool m_cond_mix_ok; + bool m_lambda_mix_ok; //! Mode indicator for transport models -- currently unused. int m_mode; diff --git a/Cantera/src/transport/LiquidTransportData.cpp b/Cantera/src/transport/LiquidTransportData.cpp index 25dcd3c82..496fb7c95 100644 --- a/Cantera/src/transport/LiquidTransportData.cpp +++ b/Cantera/src/transport/LiquidTransportData.cpp @@ -10,6 +10,7 @@ */ #include "LiquidTransportData.h" +using namespace std; namespace Cantera { @@ -74,6 +75,11 @@ namespace Cantera { speciesName = right.speciesName; hydroRadius = right.hydroRadius; viscosity = right.viscosity; + ionConductivity = right.ionConductivity; + mobilityRatio = right.mobilityRatio; + mobRatIndex = right.mobRatIndex; + selfDiffusion = right.selfDiffusion; + selfDiffIndex = right.selfDiffIndex; thermalCond = right.thermalCond; electCond = right.electCond; speciesDiffusivity = right.speciesDiffusivity; @@ -167,7 +173,9 @@ namespace Cantera { thermo_t* thermo ) : LTPspecies( propNode, name, tp_ind, thermo) { - m_model = LTR_MODEL_ARRHENIUS; + m_model = LTR_MODEL_ARRHENIUS; + m_temp = 0.0; + m_prop = 0.0; doublereal A_k, n_k, Tact_k; getArrhenius(propNode, A_k, n_k, Tact_k); @@ -268,14 +276,16 @@ namespace Cantera { thermo_t* thermo ) : LTPspecies( propNode, name, tp_ind, thermo) { - m_model = LTR_MODEL_POLY; + m_model = LTR_MODEL_POLY; + m_temp = 0.0; + m_prop = 0.0; - getFloatArray(propNode, m_coeffs, true); // if units labeled, convert Angstroms -> meters + getFloatArray(propNode, m_coeffs, "true", "toSI"); - if (m_coeffs[0] <= 0.0) { + /* if (m_coeffs[0] <= 0.0) { throw LTPError("negative or zero " + propNode.name() ); - } + }*/ } //! Copy constructor @@ -308,13 +318,89 @@ namespace Cantera { doublereal LTPspecies_Poly::getSpeciesTransProp( ) { doublereal t = m_thermo->temperature(); - if (t != m_temp) { - double tempN = 1.0; + if (t != m_temp) { + m_temp=t; + double tempN = 1.0; for (int i = 0; i < (int) m_coeffs.size() ; i++) { m_prop += m_coeffs[i] * tempN; + //cout << "m_coeff = " <temperature(); + if (t != m_temp) { + m_temp=t; + m_prop=m_coeffs[0]; + double tempN = 1.0; + for (int i = 1; i < (int) m_coeffs.size() ; i++) { + tempN *= m_temp; + m_prop *= exp(m_coeffs[i] * tempN); + //cout << "m_coeff = " < + * 0.6, -15.0e-5 + * + * + * + * + * \endverbatim + */ + class LTPspecies_ExpT : public LTPspecies{ + + public: + + LTPspecies_ExpT( const XML_Node &propNode, + std::string name, + TransportPropertyList tp_ind, + thermo_t* thermo ); + + //! Copy constructor + LTPspecies_ExpT( const LTPspecies_ExpT &right ); + + //! Assignment operator + LTPspecies_ExpT& operator=(const LTPspecies_ExpT& right ); + + virtual ~LTPspecies_ExpT( ) { } + + //! Returns the pure species tranport property + /*! + * The pure species transport property (i.e. pure species viscosity) + * is returned. Any temperature and composition dependence will be + * adjusted internally according to the information provided. + */ + doublereal getSpeciesTransProp( ); + + protected: + + //! temperature from thermo object + doublereal m_temp; + + //! most recent evaluation of transport property + doublereal m_prop; + + //! Internal model to adjust species-specific properties for composition. + /** Currently just a place holder, but this method could take + * the composition from the thermo object and adjust coefficients + * accoding to some unspecified model. + */ + void adjustCoeffsForComposition( ){ } + }; + + } #endif diff --git a/Cantera/src/transport/LiquidTransportParams.cpp b/Cantera/src/transport/LiquidTransportParams.cpp index 0caeabc59..41cebab6a 100644 --- a/Cantera/src/transport/LiquidTransportParams.cpp +++ b/Cantera/src/transport/LiquidTransportParams.cpp @@ -10,6 +10,11 @@ */ #include "LiquidTransportParams.h" +#include +#include "IonsFromNeutralVPSSTP.h" +#include "MargulesVPSSTP.h" +#include +using namespace std; namespace Cantera { @@ -49,19 +54,44 @@ namespace Cantera { { } - + LiquidTranInteraction::~LiquidTranInteraction(){ + int kmax = m_Aij.size(); + for ( int k = 0; k < kmax; k++) + if ( m_Aij[k] ) delete m_Aij[k]; + kmax = m_Bij.size(); + for ( int k = 0; k < kmax; k++) + if ( m_Bij[k] ) delete m_Bij[k]; + kmax = m_Hij.size(); + for ( int k = 0; k < kmax; k++) + if ( m_Hij[k] ) delete m_Hij[k]; + kmax = m_Sij.size(); + for ( int k = 0; k < kmax; k++) + if ( m_Sij[k] ) delete m_Sij[k]; + } void LiquidTranInteraction::init( const XML_Node &compModelNode, thermo_t* thermo ) { + + doublereal poly0; m_thermo = thermo; int nsp = thermo->nSpecies(); - m_Aij.resize( nsp, nsp, 0.0 ); m_Dij.resize( nsp, nsp, 0.0 ); m_Eij.resize( nsp, nsp, 0.0 ); - m_Sij.resize( nsp, nsp, 0.0 ); + /* + m_Aij.resize( nsp); + m_Bij.resize( nsp); + m_Hij.resize( nsp); + m_Sij.resize( nsp); + for ( int k = 0; k < nsp; k++ ){ + (*m_Aij[k]).resize( nsp, nsp, 0.0); + (*m_Bij[k]).resize( nsp, nsp, 0.0); + (*m_Hij[k]).resize( nsp, nsp, 0.0); + (*m_Sij[k]).resize( nsp, nsp, 0.0); + } + */ std::string speciesA; std::string speciesB; @@ -83,10 +113,10 @@ namespace Cantera { if ( jSpecies < 0 ) throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", "Unknown species " + speciesB ); - if ( xmlChild.hasChild( "Aij" ) ) { + /* if ( xmlChild.hasChild( "Aij" ) ) { m_Aij(iSpecies,jSpecies) = getFloat( xmlChild, "Aij", "toSI" ); m_Aij(jSpecies,iSpecies) = m_Aij(iSpecies,jSpecies) ; - } + }*/ if ( xmlChild.hasChild( "Eij" ) ) { m_Eij(iSpecies,jSpecies) = getFloat( xmlChild, "Eij", "actEnergy" ); @@ -94,11 +124,73 @@ namespace Cantera { m_Eij(jSpecies,iSpecies) = m_Eij(iSpecies,jSpecies) ; } - if ( xmlChild.hasChild( "Sij" ) ) { + if ( xmlChild.hasChild( "Aij" ) ) { + vector_fp poly; + poly0 = getFloat( poly, xmlChild, "Aij", "toSI" ); + if ( !poly.size() ) poly.push_back(poly0); + while (m_Aij.size()resize( nsp, nsp, 0.0); + m_Aij.push_back(aTemp); + } + for( int i=0; i<(int)poly.size(); i++ ){ + (*m_Aij[i])(iSpecies,jSpecies) = poly[i]; + //(*m_Aij[i])(jSpecies,iSpecies) = (*m_Aij[i])(iSpecies,jSpecies) ; + } + } + + if ( xmlChild.hasChild( "Bij" ) ) { + vector_fp poly; + poly0 = getFloat( poly, xmlChild, "Bij", "toSI" ); + if ( !poly.size() ) poly.push_back(poly0); + while (m_Bij.size()resize( nsp, nsp, 0.0); + m_Bij.push_back(bTemp); + } + for( int i=0; i<(int)poly.size(); i++ ){ + (*m_Bij[i])(iSpecies,jSpecies) = poly[i]; + //(*m_Bij[i])(jSpecies,iSpecies) = (*m_Bij[i])(iSpecies,jSpecies) ; + } + } + + if ( xmlChild.hasChild( "Hij" ) ) { + vector_fp poly; + poly0 = getFloat( poly, xmlChild, "Hij", "actEnergy" ); + if ( !poly.size() ) poly.push_back(poly0); + while (m_Hij.size()resize( nsp, nsp, 0.0); + m_Hij.push_back(hTemp); + } + for( int i=0; i<(int)poly.size(); i++ ){ + (*m_Hij[i])(iSpecies,jSpecies) = poly[i]; + (*m_Hij[i])(iSpecies,jSpecies) /= GasConstant; + //(*m_Hij[i])(jSpecies,iSpecies) = (*m_Hij[i])(iSpecies,jSpecies) ; + } + } + + if ( xmlChild.hasChild( "Sij" ) ) { + vector_fp poly; + poly0 = getFloat( poly, xmlChild, "Sij", "actEnergy" ); + if ( !poly.size() ) poly.push_back(poly0); + while (m_Sij.size()resize( nsp, nsp, 0.0); + m_Sij.push_back(sTemp); + } + for( int i=0; i<(int)poly.size(); i++ ){ + (*m_Sij[i])(iSpecies,jSpecies) = poly[i]; + (*m_Sij[i])(iSpecies,jSpecies) /= GasConstant; + //(*m_Sij[i])(jSpecies,iSpecies) = (*m_Sij[i])(iSpecies,jSpecies) ; + } + } + + /*0 if ( xmlChild.hasChild( "Sij" ) ) { m_Sij(iSpecies,jSpecies) = getFloat( xmlChild, "Sij", "toSI" ); m_Sij(iSpecies,jSpecies) /= GasConstant; - m_Sij(jSpecies,iSpecies) = m_Sij(iSpecies,jSpecies) ; - } + //m_Sij(jSpecies,iSpecies) = m_Sij(iSpecies,jSpecies) ; + }*/ if ( xmlChild.hasChild( "Dij" ) ) { m_Dij(iSpecies,jSpecies) = getFloat( xmlChild, "Dij", "toSI" ); @@ -121,7 +213,9 @@ namespace Cantera { m_thermo = right.m_thermo; //m_trParam = right.m_trParam; m_Aij = right.m_Aij; + m_Bij = right.m_Bij; m_Eij = right.m_Eij; + m_Hij = right.m_Hij; m_Sij = right.m_Sij; m_Dij = right.m_Dij; } @@ -173,6 +267,7 @@ namespace Cantera { doublereal LTI_Solvent::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); doublereal molefracs[nsp]; m_thermo->getMoleFractions(molefracs); @@ -181,23 +276,32 @@ namespace Cantera { //if weightings are specified, use those if ( speciesWeight ) { for ( int k = 0; k < nsp; k++) { - //presume that the weighting is set to 1.0 for solvent and 0.0 for everything else. - value += speciesValues[k] * speciesWeight[k]; + molefracs[k] = molefracs[k]; + // should be: molefracs[k] = molefracs[k]*speciesWeight[k]; for consistency, but weight(solvent)=1? } } else { - //This does not follow directly a solvent model + throw CanteraError("LTI_Solvent::getMixTransProp","You should be specifying the speciesWeight"); + /* //This does not follow directly a solvent model //although if the solvent mole fraction is dominant //and the other species values are given or zero, //it should work. for ( int k = 0; k < nsp; k++) { value += speciesValues[k] * molefracs[k]; - } + }*/ } - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += molefracs[i] * molefracs[j] * m_Aij(i,j) ; + for ( int i = 0; i < nsp; i++ ){ + //presume that the weighting is set to 1.0 for solvent and 0.0 for everything else. + value += speciesValues[i] * speciesWeight[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i],k); + } + for ( int k = 0; k < (int)m_Bij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i],k); + } } + } return value; } @@ -206,19 +310,29 @@ namespace Cantera { doublereal LTI_Solvent::getMixTransProp( std::vector LTPptrs ) { int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); doublereal molefracs[nsp]; m_thermo->getMoleFractions(molefracs); doublereal value = 0; for ( int k = 0; k < nsp; k++) { - //presume that the weighting is set to 1.0 for solvent and 0.0 for everything else. - value += LTPptrs[k]->getSpeciesTransProp() * LTPptrs[k]->getMixWeight( ) ; + molefracs[k] = molefracs[k]; + // should be: molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight(); for consistency, but weight(solvent)=1? } - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += molefracs[i] * molefracs[j] * m_Aij(i,j) ; + for ( int i = 0; i < nsp; i++ ){ + //presume that the weighting is set to 1.0 for solvent and 0.0 for everything else. + value += LTPptrs[i]->getSpeciesTransProp() * LTPptrs[i]->getMixWeight(); + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i],k); + } + for ( int k = 0; k < (int)m_Bij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i],k); + } + } + } return value; } @@ -229,105 +343,6 @@ namespace Cantera { doublereal LTI_MoleFracs::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { - int nsp = m_thermo->nSpecies(); - doublereal molefracs[nsp]; - m_thermo->getMoleFractions(molefracs); - - doublereal value = 0; - - //if weightings are specified, use those - if ( speciesWeight ) { - for ( int k = 0; k < nsp; k++) { - value += speciesValues[k] * speciesWeight[k] * molefracs[k]; - } - } - else { - for ( int k = 0; k < nsp; k++) { - value += speciesValues[k] * molefracs[k]; - } - } - - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += molefracs[i] * molefracs[j] * m_Aij(i,j) ; - - return value; - } - - - doublereal LTI_MoleFracs::getMixTransProp( std::vector LTPptrs ) { - - int nsp = m_thermo->nSpecies(); - doublereal molefracs[nsp]; - m_thermo->getMoleFractions(molefracs); - - doublereal value = 0; - - for ( int k = 0; k < nsp; k++) { - value += LTPptrs[k]->getSpeciesTransProp() * LTPptrs[k]->getMixWeight( ) * molefracs[k]; - } - - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += molefracs[i] * molefracs[j] * m_Aij(i,j) ; - - return value; - } - - - doublereal LTI_MassFracs::getMixTransProp(doublereal *speciesValues, doublereal *speciesWeight ) { - - int nsp = m_thermo->nSpecies(); - doublereal massfracs[nsp]; - m_thermo->getMassFractions(massfracs); - - doublereal value = 0; - - //if weightings are specified, use those - if ( speciesWeight ) { - for ( int k = 0; k < nsp; k++) { - value += speciesValues[k] * speciesWeight[k] * massfracs[k]; - } - } - else { - for ( int k = 0; k < nsp; k++) { - value += speciesValues[k] * massfracs[k]; - } - } - - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += massfracs[i] * massfracs[j] * m_Aij(i,j) ; - - return value; - } - - - doublereal LTI_MassFracs::getMixTransProp( std::vector LTPptrs ) { - - int nsp = m_thermo->nSpecies(); - doublereal massfracs[nsp]; - m_thermo->getMassFractions(massfracs); - - doublereal value = 0; - - for ( int k = 0; k < nsp; k++) { - value += LTPptrs[k]->getSpeciesTransProp() * LTPptrs[k]->getMixWeight( ) * massfracs[k]; - } - - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += massfracs[i] * massfracs[j] * m_Aij(i,j) ; - - return value; - } - - - - - - doublereal LTI_Log_MoleFracs::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { - int nsp = m_thermo->nSpecies(); doublereal temp = m_thermo->temperature(); doublereal molefracs[nsp]; @@ -338,19 +353,157 @@ namespace Cantera { //if weightings are specified, use those if ( speciesWeight ) { for ( int k = 0; k < nsp; k++) { - value += log( speciesValues[k] ) * speciesWeight[k] * molefracs[k]; + molefracs[k] = molefracs[k]*speciesWeight[k]; } } else { - for ( int k = 0; k < nsp; k++) { - value += log( speciesValues[k] ) * molefracs[k]; + throw CanteraError("LTI_MoleFracs::getMixTransProp","You should be specifying the speciesWeight"); + } + + for ( int i = 0; i < nsp; i++ ){ + value += speciesValues[i] * molefracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i],k); + } + for ( int k = 0; k < (int)m_Bij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i],k); + } } } - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += molefracs[i] * molefracs[j] - * ( m_Sij(i,j) + m_Eij(i,j) / temp ); + return value; + } + + + doublereal LTI_MoleFracs::getMixTransProp( std::vector LTPptrs ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + for ( int k = 0; k < nsp; k++) { + molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight(); + } + + for ( int i = 0; i < nsp; i++ ){ + value += LTPptrs[i]->getSpeciesTransProp() * molefracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i],k); + } + for ( int k = 0; k < (int)m_Bij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i],k); + } + } + } + return value; + } + + + doublereal LTI_MassFracs::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + doublereal massfracs[nsp]; + m_thermo->getMassFractions(massfracs); + + doublereal value = 0; + + //if weightings are specified, use those + if ( speciesWeight ) { + for ( int k = 0; k < nsp; k++) { + massfracs[k] = massfracs[k]*speciesWeight[k]; + } + } + else { + throw CanteraError("LTI_MassFracs::getMixTransProp","You should be specifying the speciesWeight"); + } + + for ( int i = 0; i < nsp; i++ ){ + value += speciesValues[i] * massfracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += massfracs[i]*massfracs[j]*(*m_Aij[k])(i,j)*pow(massfracs[i],k); + } + for ( int k = 0; k < (int)m_Bij.size(); k++ ){ + value += massfracs[i]*massfracs[j]*(*m_Bij[k])(i,j)*temp*pow(massfracs[i],k); + } + } + } + + return value; + } + + + doublereal LTI_MassFracs::getMixTransProp( std::vector LTPptrs ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + doublereal massfracs[nsp]; + m_thermo->getMassFractions(massfracs); + + doublereal value = 0; + + for ( int k = 0; k < nsp; k++) { + massfracs[k] = massfracs[k]*LTPptrs[k]->getMixWeight(); + } + + for ( int i = 0; i < nsp; i++ ){ + value += LTPptrs[i]->getSpeciesTransProp() * massfracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += massfracs[i]*massfracs[j]*(*m_Aij[k])(i,j)*pow(massfracs[i],k); + } + for ( int k = 0; k < (int)m_Bij.size(); k++ ){ + value += massfracs[i]*massfracs[j]*(*m_Bij[k])(i,j)*temp*pow(massfracs[i],k); + } + } + } + + return value; + } + + + + + doublereal LTI_Log_MoleFracs::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + + + doublereal value = 0; + + //if weightings are specified, use those + if ( speciesWeight ) { + for ( int k = 0; k < nsp; k++) { + molefracs[k] = molefracs[k]*speciesWeight[k]; + } + } + else{ + throw CanteraError("LTI_Log_MoleFracs::getMixTransProp","You probably should have a speciesWeight when you call getMixTransProp to convert ion mole fractions to molecular mole fractions"); + } + + for ( int i = 0; i < nsp; i++ ){ + value += log( speciesValues[i] ) * molefracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Hij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Hij[k])(i,j)/temp*pow(molefracs[i],k); + //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; + } + for ( int k = 0; k < (int)m_Sij.size(); k++ ){ + value -= molefracs[i]*molefracs[j]*(*m_Sij[k])(i,j)*pow(molefracs[i],k); + //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; + } + } + } value = exp( value ); return value; @@ -363,20 +516,36 @@ namespace Cantera { doublereal temp = m_thermo->temperature(); doublereal molefracs[nsp]; m_thermo->getMoleFractions(molefracs); + doublereal value = 0; + //if weightings are specified, use those + for ( int k = 0; k < nsp; k++) { - value += log( LTPptrs[k]->getSpeciesTransProp() ) * LTPptrs[k]->getMixWeight() * molefracs[k]; + molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight( ); } - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - value += molefracs[i] * molefracs[j] - * ( m_Sij(i,j) + m_Eij(i,j) / temp ); + for ( int i = 0; i < nsp; i++ ){ + value += log( LTPptrs[i]->getSpeciesTransProp() ) * molefracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Hij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Hij[k])(i,j)/temp*pow(molefracs[i],k); + //cout << "1 = " << molefracs[i]+molefracs[j] << endl; + //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; + } + for ( int k = 0; k < (int)m_Sij.size(); k++ ){ + value -= molefracs[i]*molefracs[j]*(*m_Sij[k])(i,j)*pow(molefracs[i],k); + //cout << "1 = " << molefracs[i]+molefracs[j] << endl; + //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; + } + } + } value = exp( value ); - + // cout << ", viscSpeciesA = " << LTPptrs[0]->getSpeciesTransProp() << endl; + //cout << ", viscSpeciesB = " << LTPptrs[1]->getSpeciesTransProp() << endl; + //cout << "value = " << value << " FINAL" << endl; return value; } @@ -422,7 +591,7 @@ namespace Cantera { return value; } - void LTI_Pairwise_Interaction::getMatrixTransProp(DenseMatrix &mat, doublereal *speciesValues) { + void LTI_Pairwise_Interaction::getMatrixTransProp( DenseMatrix &mat, doublereal *speciesValues ) { int nsp = m_thermo->nSpecies(); doublereal temp = m_thermo->temperature(); @@ -439,6 +608,205 @@ namespace Cantera { mat(i,i) = m_diagonals[i]->getSpeciesTransProp() ; } + + void LTI_StefanMaxwell_PPN::setParameters( LiquidTransportParams& trParam ) { + int nsp = m_thermo->nSpecies(); + int nBinInt = nsp*(nsp-1)/2; + //vectornSpecies(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + throw LTPmodelError( "Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense." ); + + return value; + } + + + doublereal LTI_StefanMaxwell_PPN::getMixTransProp( std::vector LTPptrs ) { + + int nsp = m_thermo->nSpecies(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + throw LTPmodelError( "Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense." ); + + return value; + } + + void LTI_StefanMaxwell_PPN::getMatrixTransProp( DenseMatrix &mat, doublereal *speciesValues ) { + //CAL + + IonsFromNeutralVPSSTP * ions_thermo = dynamic_cast(m_thermo); + int i, j, k; + int nsp = m_thermo->nSpecies(); + if (nsp != 3) throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Function may only be called with a 3-ion system"); + int nBinInt = nsp*(nsp-1)/2; + doublereal temp = m_thermo->temperature(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions( molefracs ); + vector_fp neut_molefracs; + ions_thermo->getNeutralMolecMoleFractions(neut_molefracs); + vector cation; + vector anion; + ions_thermo->getCationList(cation); + ions_thermo->getAnionList(anion); + vector speciesNames; + ions_thermo->getSpeciesNames(speciesNames); + + // Reaction Coeffs and Charges + std::vector viS(6); + std::vector charges(3); + ions_thermo->getDissociationCoeffs(viS,charges); + + if ((int)anion.size() != 1) + throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Must have one anion only for StefanMaxwell_PPN"); + if ((int)cation.size() != 2) + throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Must have two cations of equal charge for StefanMaxwell_PPN"); + if (charges[cation[0]] != charges[cation[1]]) + throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Cations must be of equal charge for StefanMaxwell_PPN"); + + /* + cout << "cation 0: " << speciesNames[cation[0]] << endl; + cout << "cation 1: " << speciesNames[cation[1]] << endl; + cout << "anion 0: " << speciesNames[anion[0]] << endl; + */ + + doublereal kappa = m_ionCondMixModel->getMixTransProp(m_ionCondSpecies); + + MargulesVPSSTP * marg_thermo = dynamic_cast (ions_thermo->neutralMoleculePhase_); + doublereal vol = m_thermo->molarVolume(); + + typedef std::vector intVec; + std::vector mobRatIndexMap; + mobRatIndexMap.resize(nsp); + for ( k = 0; k < nsp; k++ ) + mobRatIndexMap[k].resize(nsp,0); + + for ( k = 0; k < nBinInt; k++ ) { + bool missingIndex = 1; + for ( i = 0; i < nsp; i++ ) { + for ( j = 0; j < i; j++ ) { + if ( (speciesNames[i] + ":" + speciesNames[j]) == m_mobRatIndex[k] ) { + mobRatIndexMap[i][j] = k; + m_mobRatMix(i,j) = m_mobRatMixModel[k]->getMixTransProp( m_mobRatSpecies[k] ); + if ( m_mobRatMix(i,j) > 0 ) m_mobRatMix(j,i) = 1.0/m_mobRatMix(i,j); + missingIndex = 0; + break; + } + else if ( (speciesNames[j] + ":" + speciesNames[i]) == m_mobRatIndex[k] ) { + m_mobRatMix(j,i) = m_mobRatMixModel[k]->getMixTransProp( m_mobRatSpecies[k] ); + if ( m_mobRatMix(j,i) > 0 ) m_mobRatMix(i,j) = 1.0/m_mobRatMix(j,i); + missingIndex = 0; + break; + } + } + } + if ( missingIndex ) throw CanteraError("LTI_StefanMaxwell_PPN::getMixTransProp","Incorrect names for mobility ratio of " + m_mobRatIndex[k] + " rather than i.e. " + speciesNames[0] + ":" + speciesNames[1]); + } + + + for ( k = 0; k < nsp; k++ ){ + j = 0; + while (m_selfDiffIndex[k] != speciesNames[j]) { + j++; + if (j == nsp) throw CanteraError("LTI_StefanMaxwell_PPN::getMixTransProp","Incorrect names for self diffusion of " + m_selfDiffIndex[k] + " rather than i.e. " + speciesNames[0]); + } + m_selfDiffMix[j] = m_selfDiffMixModel[k]->getMixTransProp( m_selfDiffSpecies[k] ); + } + + /* + for ( i = 0; i < nsp; i++ ) { + cout << "D" << i << "* = " << m_selfDiffMix[i] << endl; + for ( j = 0; j < nsp; j++ ) { + cout << "ratio" << i << j << " = " << m_mobRatMix(i,j) << endl; + } + } + */ + + int vP = max(viS[cation[0]],viS[cation[1]]); + int vM = viS[anion[0]]; + int zP = charges[cation[0]]; + int zM = charges[anion[0]]; + doublereal xA, xB, eps; + doublereal inv_vP_vM_MutualDiff; + vector_fp dlnActCoeffdlnN; + dlnActCoeffdlnN.resize(neut_molefracs.size(),0.0); + + std::string cationIndex (4,'0'); + for ( i = 0; i < 2; i++ ) + for ( j = 0; j < 2; j++ ) + if ( viS[i*nsp+cation[j]] > 0 ) + cationIndex[i*2+j] = '1'; + + if ( (cationIndex == "1001") | (cationIndex == "0110") ) { + xA = neut_molefracs[cation[0]]; + xB = neut_molefracs[cation[1]]; + eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0])); + marg_thermo->getdlnActCoeffdlnN(&dlnActCoeffdlnN[0]); + inv_vP_vM_MutualDiff = (xA*(1+dlnActCoeffdlnN[cation[1]])/m_selfDiffMix[cation[1]]+xB*(1+dlnActCoeffdlnN[cation[0]])/m_selfDiffMix[cation[0]]); + //marg_thermo->getdlnActCoeffdlnX(&dlnActCoeffdlnN[0]); + //inv_vP_vM_MutualDiff = (xA*(1+dlnActCoeffdlnN[cation[1]])/m_selfDiffMix[cation[1]]+xB*(1+dlnActCoeffdlnN[cation[0]])/m_selfDiffMix[cation[0]]); + } + else + throw CanteraError("LTI_StefanMaxwell_PPN::getMixTransProp","Dissociation reactions don't make sense: cationIndex = " + cationIndex); + + mat.resize( nsp, nsp, 0.0 ); + mat(cation[0],cation[1]) = mat(cation[1],cation[0]) = (1+vM/vP)*(1+eps*xB)*(1-eps*xA)*inv_vP_vM_MutualDiff-zP*zP*Faraday*Faraday/GasConstant/temp/kappa/vol; + mat(cation[0],anion[0]) = mat(anion[0],cation[0]) = (1+vP/vM)*(-eps*xB*(1-eps*xA)*inv_vP_vM_MutualDiff)-zP*zM*Faraday*Faraday/GasConstant/temp/kappa/vol; +mat(cation[1],anion[0]) = mat(anion[0],cation[1]) = (1+vP/vM)*(eps*xA*(1+eps*xB)*inv_vP_vM_MutualDiff)-zP*zM*Faraday*Faraday/GasConstant/temp/kappa/vol; + + + for ( i = 0; i < nsp; i++ ) { + for ( j = 0; j < nsp; j++ ) { + mat(i,j) = 1.0/mat(i,j); + //cout << "D" << i << j << " = " << mat(i,j) << endl; + } + } + } + doublereal LTI_StokesEinstein::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { @@ -503,6 +871,61 @@ namespace Cantera { delete viscSpec; } + doublereal LTI_MoleFracs_ExpT::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + //if weightings are specified, use those + if ( speciesWeight ) { + for ( int k = 0; k < nsp; k++) { + molefracs[k] = molefracs[k]*speciesWeight[k]; + } + } + else { + throw CanteraError("LTI_MoleFracs_ExpT::getMixTransProp","You should be specifying the speciesWeight"); + } + + for ( int i = 0; i < nsp; i++ ){ + value += speciesValues[i] * molefracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i],k)*exp((*m_Bij[k])(i,j)*temp); + } + } + } + + return value; + } + + + doublereal LTI_MoleFracs_ExpT::getMixTransProp( std::vector LTPptrs ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + for ( int k = 0; k < nsp; k++) { + molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight(); + } + + for ( int i = 0; i < nsp; i++ ){ + value += LTPptrs[i]->getSpeciesTransProp() * molefracs[i]; + for ( int j = 0; j < nsp; j++ ){ + for ( int k = 0; k < (int)m_Aij.size(); k++ ){ + value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i],k)*exp((*m_Bij[k])(i,j)*temp); + } + } + } + return value; + } diff --git a/Cantera/src/transport/LiquidTransportParams.h b/Cantera/src/transport/LiquidTransportParams.h index 0757978a7..cf9570d97 100644 --- a/Cantera/src/transport/LiquidTransportParams.h +++ b/Cantera/src/transport/LiquidTransportParams.h @@ -33,57 +33,62 @@ namespace Cantera { }; - //! Composition dependence type for liquid mixture transport properties - /*! - * Types of temperature dependencies: - * - 0 - Mixture calculations with this property are not allowed - * - 1 - Use solvent (species 0) properties - * - 2 - Properties weighted linearly by mole fractions - * - 3 - Properties weighted linearly by mass fractions - * - 4 - Properties weighted logarithmically by mole fractions (interaction energy weighting) - * - 5 - Interactions given pairwise between each possible species (i.e. D_ij) - * - * \verbatim - * - * - * - * - * LiCl(L) - * KCl(L) - * -1.0 - * 1.0E-1 - * - * - * - * - * - * - * Li+ - * K+ - * 1.5 - * - * - * K+ - * Cl- - * 1.0 - * - * - * Li+ - * Cl- - * 1.2 - * - * - * - * - * - * - * - * - * - * - * \endverbatim - * - */ + //! Composition dependence type for liquid mixture transport properties + /*! + * Types of temperature dependencies: + * - 0 - Mixture calculations with this property are not allowed + * - 1 - Use solvent (species 0) properties + * - 2 - Properties weighted linearly by mole fractions + * - 3 - Properties weighted linearly by mass fractions + * - 4 - Properties weighted logarithmically by mole fractions (interaction energy weighting) + * - 5 - Interactions given pairwise between each possible species (i.e. D_ij) + * + * \verbatim + * + * + * + * + * LiCl(L) + * KCl(L) + * -1.0 + * 1.0E-1 + * -or- + * 1.0E-1, 0.001 0.01 + * + * -same form for Hij,Aij,Bij- + * + * + * + * + * + * + * Li+ + * K+ + * 1.5 + * + * + * K+ + * Cl- + * 1.0 + * + * + * Li+ + * Cl- + * 1.2 + * + * + * + * + * + * + * + * + * + * + * \endverbatim + * + */ + enum LiquidTranMixingModel { LTI_MODEL_NOTSET=-1, LTI_MODEL_NONE, @@ -92,7 +97,9 @@ namespace Cantera { LTI_MODEL_MASSFRACS, LTI_MODEL_LOG_MOLEFRACS, LTI_MODEL_PAIRWISE_INTERACTION, - LTI_MODEL_STOKES_EINSTEIN + LTI_MODEL_STEFANMAXWELL_PPN, + LTI_MODEL_STOKES_EINSTEIN, + LTI_MODEL_MOLEFRACS_EXPT }; @@ -132,7 +139,7 @@ namespace Cantera { LiquidTranInteraction& operator=( const LiquidTranInteraction &right ); //! destructor - virtual ~LiquidTranInteraction() { ; } + virtual ~LiquidTranInteraction(); //! initialize LiquidTranInteraction objects with thermo and XML node /** @@ -172,14 +179,21 @@ namespace Cantera { //LiquidTransportParams* m_trParam; - //! Matrix of interactions (no temperature dependence, dimensionless) - DenseMatrix m_Aij; + //! Matrix of interaction coefficients for polynomial in molefraction*weight of speciesA (no temperature dependence, dimensionless) + std::vector m_Aij; + + //! Matrix of interaction coefficients for polynomial in molefraction*weight of speciesA (linear temperature dependence, units 1/K) + std::vector m_Bij; //! Matrix of interactions (in energy units, 1/RT temperature dependence) DenseMatrix m_Eij; - //! Matrix of interactions (in entropy units, divided by R) - DenseMatrix m_Sij; + //! Matrix of interaction coefficients for polynomial in molefraction*weight of speciesA (in energy units, 1/RT temperature dependence) + std::vector m_Hij; + + //! Matrix of interaction coefficients for polynomial in molefraction*weight of speciesA (in entropy units, divided by R) + std::vector m_Sij; + //DenseMatrix m_Sij; //! Matrix of interactions DenseMatrix m_Dij; @@ -205,6 +219,11 @@ namespace Cantera { std::vector LTData; LiquidTranInteraction* viscosity; + LiquidTranInteraction* ionConductivity; + std::vector mobilityRatio; + std::vector mobRatIndex; + std::vector selfDiffusion; + std::vector selfDiffIndex; LiquidTranInteraction* thermalCond; LiquidTranInteraction* speciesDiffusivity; LiquidTranInteraction* electCond; @@ -213,7 +232,19 @@ namespace Cantera { //! Model for species interaction effects for viscosity //! Takes enum LiquidTranMixingModel LiquidTranMixingModel model_viscosity; + + //! Model for species interaction effects for ionic conductivity + //! Takes enum LiquidTranMixingModel + LiquidTranMixingModel model_ionConductivity; + //! Model for species interaction effects for mobility ratio + //! Takes enum LiquidTranMixingModel + std::vector model_mobilityRatio; + + //! Model for species interaction effects for mobility ratio + //! Takes enum LiquidTranMixingModel + std::vector model_selfDiffusion; + //! Interaction associated with linear weighting of //! thermal conductivity. /** @@ -233,7 +264,7 @@ namespace Cantera { //! Interaction associated with linear weighting of //! thermal conductivity. /** - * This is used for either LTI_MODEL_PAIRWISE_INTERACTION. + * This is used for either LTI_MODEL_PAIRWISE_INTERACTION or LTI_MODEL_STEFANMAXWELL_PPN. * These provide species interaction coefficients associated with * the Stefan-Maxwell formulation. */ @@ -256,7 +287,8 @@ namespace Cantera { class LTI_MassFracs; class LTI_Log_MoleFracs; class LTI_Pairwise_Interaction; - + class LTI_StefanMaxwell_PPN; + class LTI_MoleFracs_ExpT; class LTI_Solvent : public LiquidTranInteraction { public: @@ -283,7 +315,7 @@ namespace Cantera { doublereal getMixTransProp( doublereal *valueSpecies, doublereal *weightSpecies = 0 ); doublereal getMixTransProp( std::vector LTPptrs ) ; - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = m_Aij; } + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } protected: @@ -325,7 +357,7 @@ namespace Cantera { doublereal getMixTransProp( doublereal *valueSpecies, doublereal *weightSpecies = 0 ); doublereal getMixTransProp( std::vector LTPptrs ) ; - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = m_Aij; } + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } protected: @@ -368,7 +400,7 @@ namespace Cantera { doublereal getMixTransProp( doublereal *valueSpecies, doublereal *weightSpecies = 0 ); doublereal getMixTransProp( std::vector LTPptrs ) ; - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = m_Aij; } + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } protected: @@ -502,13 +534,84 @@ namespace Cantera { doublereal getMixTransProp( std::vector LTPptrs ) ; void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) ; - protected: std::vector m_diagonals; }; + //! Stefan Maxwell Diffusion Coefficients can be solved for given + //! ion conductivity, mobility ratios, and self diffusion coeffs. + //! This method is only valid for a common anion mixture of two + //! salts with cations of equal charge. + /** + * This class requres you specify + * 1 - ion conductivity + * 2 - mobility ratio of the two cations + * 3 - mutual diffusion coefficient (can be approximated using + * the self diffusion coefficients of the cations + * + * Sample input for this method is + * \verbatim + * + * + * + * + * + * + * \endverbatim + * + */ + class LTI_StefanMaxwell_PPN : public LiquidTranInteraction { + + public: + LTI_StefanMaxwell_PPN( TransportPropertyList tp_ind = TP_UNKNOWN ) : + LiquidTranInteraction( tp_ind ) + { + m_model = LTI_MODEL_STEFANMAXWELL_PPN; + } + + + //! Copy constructor + // LTI_StefanMaxwell_PPN( const LTI_StefanMaxwell_PPN &right ); + + //! Assignment operator + // LTI_StefanMaxwell_PPN& operator=( const LTI_StefanMaxwell_PPN &right ); + + virtual ~LTI_StefanMaxwell_PPN( ) { } + + void setParameters( LiquidTransportParams& trParam ) ; + + //! Return the mixture transport property value. + /** + * Takes the separate species transport properties + * as input (this method does not know what + * transport property it is at this point. + */ + doublereal getMixTransProp( doublereal *valueSpecies, doublereal *weightSpecies = 0 ); + doublereal getMixTransProp( std::vector LTPptrs ) ; + + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) ; + //CAL void getMatrixTransProp( DenseMatrix &mat, LiquidTransport* lt, doublereal* speciesValues = 0 ) ; + + protected: + + doublereal m_ionCondMix; + LiquidTranInteraction * m_ionCondMixModel; + std::vector m_ionCondSpecies; + typedef std::vector LTPvector; + DenseMatrix m_mobRatMix; + std::vector m_mobRatMixModel; + std::vector m_mobRatSpecies; + std::vector m_mobRatIndex; + + std::vector m_selfDiffMixModel; + vector_fp m_selfDiffMix; + std::vector m_selfDiffSpecies; + std::vector m_selfDiffIndex; + }; + + class LTI_StokesEinstein : public LiquidTranInteraction { public: @@ -539,7 +642,6 @@ namespace Cantera { doublereal getMixTransProp( std::vector LTPptrs ) ; void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) ; - protected: std::vector m_viscosity; @@ -547,6 +649,48 @@ namespace Cantera { }; + //! Simple mole fraction weighting of transport properties + /** + * This model weights the transport property by the mole + * fractions. + * The overall formula for the mixture viscosity is + * + * \f[ \eta_{mix} = \sum_i X_i \eta_i + * + \sum_i \sum_j X_i X_j A_{i,j} \f]. + */ + class LTI_MoleFracs_ExpT : public LiquidTranInteraction { + + public: + LTI_MoleFracs_ExpT( TransportPropertyList tp_ind = TP_UNKNOWN ) : + LiquidTranInteraction( tp_ind ) + { + m_model = LTI_MODEL_MOLEFRACS_EXPT; + } + + + //! Copy constructor + // LTI_MoleFracs_ExpT( const LTI_MoleFracs_ExpT &right ); + + //! Assignment operator + // LTI_MoleFracs_ExpT& operator=( const LTI_MoleFracs_ExpT &right ); + + virtual ~LTI_MoleFracs_ExpT( ) { } + + //! Return the mixture transport property value. + /** + * Takes the separate species transport properties + * as input (this method does not know what + * transport property it is at this point. + */ + doublereal getMixTransProp( doublereal *valueSpecies, doublereal *weightSpecies = 0 ); + doublereal getMixTransProp( std::vector LTPptrs ) ; + + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } + //CAL void getMatrixTransProp( DenseMatrix &mat, LiquidTransport* lt, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } + + protected: + + }; } diff --git a/Cantera/src/transport/TransportBase.h b/Cantera/src/transport/TransportBase.h index 24ad9d137..7dd194f8c 100644 --- a/Cantera/src/transport/TransportBase.h +++ b/Cantera/src/transport/TransportBase.h @@ -28,7 +28,7 @@ #define CT_TRANSPORTBASE_H #include "ThermoPhase.h" - +#include "DenseMatrix.h" namespace Cantera { @@ -80,6 +80,8 @@ namespace Cantera { * - VB_MOLEAVG Diffusion velocities are based on the mole averaged velocities * - VB_SPECIES_0 Diffusion velocities are based on the relative motion wrt species 0 * - VB_SPECIES_1 Diffusion velocities are based on the relative motion wrt species 1 + * - VB_SPECIES_2 Diffusion velocities are based on the relative motion wrt species 2 + * - VB_SPECIES_3 Diffusion velocities are based on the relative motion wrt species 3 * * @ingroup tranprops */ @@ -97,6 +99,10 @@ namespace Cantera { const VelocityBasis VB_SPECIES_0 = 0; //! Diffusion velocities are based on the relative motion wrt species 1 const VelocityBasis VB_SPECIES_1 = 1; + //! Diffusion velocities are based on the relative motion wrt species 2 + const VelocityBasis VB_SPECIES_2 = 2; + //! Diffusion velocities are based on the relative motion wrt species 3 + const VelocityBasis VB_SPECIES_3 = 3; //@} //! Base class for transport property managers. @@ -308,6 +314,68 @@ namespace Cantera { virtual doublereal bulkViscosity() { return err("bulkViscosity"); } + /** + * The ionic conducitivity in 1/ohm/m. + */ + virtual doublereal ionConductivity() + { return err("ionConductivity"); } + + //! Returns the pure species ionic conducitivity + /*! + * The units are 1/ohm/m and the length is the number of species + * + * @param ionCond Vector of ionic conductivities + */ + virtual void getSpecoesIonConductivity(doublereal* const ionCond) + { err("getSpeciesIonConductivity"); } + + + + /** + * The mobility ratio. + */ + virtual void mobilityRatio(vector_fp& mobRat, std::vector& mobRatIndex) + { err("mobilityRatio"); } + virtual void mobilityRatio(double* mobRat, std::vector& mobRatIndex) + { err("mobilityRatio"); } + + //! Returns the pure species limit of the mobility ratios + /*! + * The value is dimensionless and the length is the number of species + * + * @param mobRat Vector of mobility ratios + */ + virtual void getSpeciesMobilityRatio(DenseMatrix& mobRat, std::vector& mobRatIndex) + { err("getSpeciesMobilityRatio"); } + virtual void getSpeciesMobilityRatio(double** mobRat, std::vector& mobRatIndex) + { err("getSpeciesMobilityRatio"); } + + //! Returns the self diffusion coefficients in the solution + /*! + * The self diffusion calculation is handled by subclasses of + * LiquidTranInteraction as specified in the input file. + * These in turn employ subclasses of LTPspecies to + * determine the individual species self diffusion coeffs. + */ + virtual void selfDiffusion(vector_fp& selfDiff, std::vector& selfDiffIndex) + { err("selfDiffusion"); } + virtual void selfDiffusion(double* selfDiff, std::vector& selfDiffIndex) + { err("selfdiffusion"); } + + //! Returns the pure species self diffusion in solution of each species + /*! + * The pure species molar volumes are evaluated using the + * appropriate subclasses of LTPspecies as specified in the + * input file. + * + * @param selfDiff array of length "number of species" + * to hold returned self diffusion coeffs. + */ + virtual void getSpeciesSelfDiffusion(DenseMatrix& selfDiff, std::vector& selfDiffIndex) + { err("getSpeciesSelfDiffusion"); } + virtual void getSpeciesSelfDiffusion(double** selfDiff, std::vector& selfDiffIndex) + { err("getSpeciesSelfDiffusion"); } + //! Returns the mixture thermal conductivity in W/m/K. /*! * Units are in W / m K or equivalently kg m / s3 K diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index a69ba0d77..f000345e2 100644 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -232,6 +232,9 @@ namespace Cantera { //m_models["Radiative"] = cRadiative; m_tranPropMap["viscostiy"] = TP_VISCOSITY; + m_tranPropMap["ionConductivity"] = TP_IONCONDUCTIVITY; + m_tranPropMap["mobilityRatio"] = TP_MOBILITYRATIO; + m_tranPropMap["selfDiffusion"] = TP_SELFDIFFUSION; m_tranPropMap["thermalConductivity"] = TP_THERMALCOND; m_tranPropMap["speciesDiffusivity"] = TP_DIFFUSIVITY; m_tranPropMap["hydrodynamicRadius"] = TP_HYDRORADIUS; @@ -241,6 +244,7 @@ namespace Cantera { m_LTRmodelMap["constant"] = LTR_MODEL_CONSTANT; m_LTRmodelMap["arrhenius"] = LTR_MODEL_ARRHENIUS; m_LTRmodelMap["coeffs"] = LTR_MODEL_POLY; + m_LTRmodelMap["exptemp"] = LTR_MODEL_EXPT; m_LTImodelMap[""] = LTI_MODEL_NOTSET; m_LTImodelMap["none"] = LTI_MODEL_NONE; @@ -249,6 +253,8 @@ namespace Cantera { m_LTImodelMap["massFractions"] = LTI_MODEL_MASSFRACS; m_LTImodelMap["logMoleFractions"] = LTI_MODEL_LOG_MOLEFRACS; m_LTImodelMap["pairwiseInteraction"] = LTI_MODEL_PAIRWISE_INTERACTION; + m_LTImodelMap["stefanMaxwell_PPN"] = LTI_MODEL_STEFANMAXWELL_PPN; + m_LTImodelMap["moleFractionsExpT"] = LTI_MODEL_MOLEFRACS_EXPT; } /** @@ -312,6 +318,12 @@ namespace Cantera { tp_ind, thermo ); break; + case LTR_MODEL_EXPT: + ltps = new LTPspecies_ExpT( trNode, + name, + tp_ind, + thermo ); + break; default: throw CanteraError("newLTP","unknown transport model: " + model ); ltps = new LTPspecies( trNode, @@ -359,11 +371,20 @@ namespace Cantera { lti->init( trNode, thermo ); lti->setParameters( trParam ); break; + case LTI_MODEL_STEFANMAXWELL_PPN: + lti = new LTI_StefanMaxwell_PPN( tp_ind ); + lti->init( trNode, thermo ); + lti->setParameters( trParam ); + break; case LTI_MODEL_STOKES_EINSTEIN: lti = new LTI_StokesEinstein( tp_ind ); lti->init( trNode, thermo ); lti->setParameters( trParam ); break; + case LTI_MODEL_MOLEFRACS_EXPT: + lti = new LTI_MoleFracs_ExpT( tp_ind ); + lti->init( trNode, thermo ); + break; default: // throw CanteraError("newLTI","unknown transport model: " + model ); lti = new LiquidTranInteraction( tp_ind ); @@ -966,7 +987,7 @@ namespace Cantera { for (int iChild = 0; iChild < num; iChild++) { XML_Node &xmlChild = trNode.child(iChild); std::string nodeName = xmlChild.name(); - + switch ( m_tranPropMap[nodeName] ) { case TP_VISCOSITY: data.viscosity = newLTP( xmlChild, @@ -974,23 +995,63 @@ namespace Cantera { m_tranPropMap[nodeName], trParam.thermo ); break; + case TP_IONCONDUCTIVITY: + data.ionConductivity = newLTP( xmlChild, + name, + m_tranPropMap[nodeName], + trParam.thermo ); + break; + case TP_MOBILITYRATIO: + { + int iSpec; + int numSpec = xmlChild.nChildren(); + data.mobRatIndex.resize(numSpec); + data.mobilityRatio.resize(numSpec); + for (iSpec = 0; iSpec< numSpec; iSpec++){ + XML_Node &propSpecNode = xmlChild.child(iSpec); + std::string specName = propSpecNode.name(); + data.mobRatIndex[iSpec] = specName; + data.mobilityRatio[iSpec] = newLTP( propSpecNode, + name, + m_tranPropMap[nodeName], + trParam.thermo ); + }; + }; + break; + case TP_SELFDIFFUSION: + { + int iSpec; + int numSpec = xmlChild.nChildren(); + data.selfDiffIndex.resize(numSpec); + data.selfDiffusion.resize(numSpec); + for (iSpec = 0; iSpec< numSpec; iSpec++){ + XML_Node &propSpecNode = xmlChild.child(iSpec); + std::string specName = propSpecNode.name(); + data.selfDiffIndex[iSpec] = specName; + data.selfDiffusion[iSpec] = newLTP( propSpecNode, + name, + m_tranPropMap[nodeName], + trParam.thermo ); + }; + }; + break; case TP_THERMALCOND: data.thermalCond = newLTP( xmlChild, - name, - m_tranPropMap[nodeName], - trParam.thermo ); + name, + m_tranPropMap[nodeName], + trParam.thermo ); break; case TP_DIFFUSIVITY: data.speciesDiffusivity = newLTP( xmlChild, - name, - m_tranPropMap[nodeName], - trParam.thermo ); + name, + m_tranPropMap[nodeName], + trParam.thermo ); break; case TP_HYDRORADIUS: data.hydroRadius = newLTP( xmlChild, - name, - m_tranPropMap[nodeName], - trParam.thermo ); + name, + m_tranPropMap[nodeName], + trParam.thermo ); break; case TP_ELECTCOND: data.electCond = newLTP( xmlChild, @@ -1005,8 +1066,8 @@ namespace Cantera { } } - - datatable[name] = data; + datatable.insert(pair(name,data)); + // datatable[name] = data; } } catch(CanteraError) { @@ -1059,12 +1120,10 @@ namespace Cantera { for (int iChild = 0; iChild < num; iChild++) { //tranTypeNode is a type of transport property like viscosity XML_Node &tranTypeNode = transportNode.child(iChild); - if (tranTypeNode.hasChild("compositionDependence")) { + std::string nodeName = tranTypeNode.name(); + if ( tranTypeNode.hasChild("compositionDependence")) { //compDepNode contains the interaction model XML_Node &compDepNode = tranTypeNode.child("compositionDependence"); - - std::string nodeName = tranTypeNode.name(); - switch (m_tranPropMap[nodeName]) { break; case TP_VISCOSITY: @@ -1072,6 +1131,43 @@ namespace Cantera { m_tranPropMap[nodeName], trParam ); break; + case TP_IONCONDUCTIVITY: + trParam.ionConductivity = newLTI( compDepNode, + m_tranPropMap[nodeName], + trParam ); + break; + case TP_MOBILITYRATIO: + { + int iSpec; + int numSpec = compDepNode.nChildren(); + trParam.mobRatIndex.resize(numSpec); + trParam.mobilityRatio.resize(numSpec); + for (iSpec = 0; iSpec< numSpec; iSpec++){ + XML_Node &propSpecNode = compDepNode.child(iSpec); + std::string specName = propSpecNode.name(); + trParam.mobRatIndex[iSpec] = specName; + trParam.mobilityRatio[iSpec] = newLTI( propSpecNode, + m_tranPropMap[nodeName], + trParam ); + }; + }; + break; + case TP_SELFDIFFUSION: + { + int iSpec; + int numSpec = compDepNode.nChildren(); + trParam.selfDiffIndex.resize(numSpec); + trParam.selfDiffusion.resize(numSpec); + for (iSpec = 0; iSpec< numSpec; iSpec++){ + XML_Node &propSpecNode = compDepNode.child(iSpec); + std::string specName = propSpecNode.name(); + trParam.selfDiffIndex[iSpec] = specName; + trParam.selfDiffusion[iSpec] = newLTI( propSpecNode, + m_tranPropMap[nodeName], + trParam ); + }; + }; + break; case TP_THERMALCOND: trParam.thermalCond = newLTI( compDepNode, m_tranPropMap[nodeName], @@ -1122,15 +1218,17 @@ namespace Cantera { trParam.velocityBasis_ = trParam.thermo->speciesIndex(velocityBasis) ; else { int linenum; - throw TransportDBError( linenum, "Unknown attribute " + velocityBasis + " for node. "); + throw TransportDBError( linenum, "Unknown attribute \"" + velocityBasis + "\" for node. "); } - } } } - catch(CanteraError) { - ; + catch (CanteraError) { + showErrors(std::cout); } + //catch(CanteraError) { + // ; + //} return; }