From 497e69d0d0d69657502cbbad1d895fad4cf5320d Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 10 Apr 2010 02:35:18 +0000 Subject: [PATCH] Worked on doxygen updates Split a file in two in order to make the directory more readable. Started working on standardizing the files. --- .../src/transport/LiquidTranInteraction.cpp | 864 +++++++++++++++++ Cantera/src/transport/LiquidTranInteraction.h | 685 +++++++++++++ .../src/transport/LiquidTransportParams.cpp | 902 ++---------------- Cantera/src/transport/LiquidTransportParams.h | 645 +------------ Cantera/src/transport/Makefile.in | 4 +- Cantera/src/transport/TransportFactory.cpp | 97 +- Cantera/src/transport/TransportFactory.h | 27 +- Cantera/src/transport/TransportParams.h | 7 + 8 files changed, 1724 insertions(+), 1507 deletions(-) create mode 100644 Cantera/src/transport/LiquidTranInteraction.cpp create mode 100644 Cantera/src/transport/LiquidTranInteraction.h diff --git a/Cantera/src/transport/LiquidTranInteraction.cpp b/Cantera/src/transport/LiquidTranInteraction.cpp new file mode 100644 index 000000000..854b7a7ac --- /dev/null +++ b/Cantera/src/transport/LiquidTranInteraction.cpp @@ -0,0 +1,864 @@ +/** + * @file LiquidTransportParams.cpp + * Source code for liquid mixture transport property evaluations. + */ +/* + * Latest Checkin: + * $Author: hkmoffa $ + * $Date: 2010-04-02 11:40:52 -0600 (Fri, 02 Apr 2010) $ + * $Revision: 429 $ + */ + +#include "LiquidTransportParams.h" +#include +#include "IonsFromNeutralVPSSTP.h" +#include "MargulesVPSSTP.h" +#include +using namespace std; + + +namespace Cantera { + + /** + * Exception thrown if an error is encountered while reading the + * transport database. + */ + class LTPError : public CanteraError { + public: + LTPError(std::string msg ) + : CanteraError("LTPspecies", + "error parsing transport data: " + + msg + "\n") {} + }; + + /** + * Exception thrown if an error is encountered while reading the + * transport database. + */ + class LTPmodelError : public CanteraError { + public: + LTPmodelError(std::string msg ) + : CanteraError("LTPspecies", + "error parsing transport data: " + + msg + "\n") {} + }; + + + // Constructor + /* + * @param tp_ind Index indicating transport property type (i.e. viscosity) + */ + LiquidTranInteraction::LiquidTranInteraction(TransportPropertyList tp_ind ) : + m_model(LTI_MODEL_NOTSET), + m_property(tp_ind) + { + } + + 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_Dij.resize(nsp, nsp, 0.0 ); + m_Eij.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; + + int num = compModelNode.nChildren(); + for (int iChild = 0; iChild < num; iChild++) { + XML_Node &xmlChild = compModelNode.child(iChild); + std::string nodeName = lowercase(xmlChild.name() ); + if (nodeName != "interaction" ) + throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", + "expected element and got <" + nodeName + ">" ); + speciesA = xmlChild.attrib("speciesA"); + speciesB = xmlChild.attrib("speciesB"); + int iSpecies = m_thermo->speciesIndex(speciesA ); + if (iSpecies < 0 ) + throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", + "Unknown species " + speciesA ); + int jSpecies = m_thermo->speciesIndex(speciesB ); + if (jSpecies < 0 ) + throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", + "Unknown species " + speciesB ); + /* 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" ); + m_Eij(iSpecies,jSpecies) /= GasConstant; + m_Eij(jSpecies,iSpecies) = m_Eij(iSpecies,jSpecies) ; + } + + 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) ; + }*/ + + if (xmlChild.hasChild("Dij" ) ) { + m_Dij(iSpecies,jSpecies) = getFloat(xmlChild, "Dij", "toSI" ); + m_Dij(jSpecies,iSpecies) = m_Dij(iSpecies,jSpecies) ; + } + } + } + + // Copy constructor + LiquidTranInteraction::LiquidTranInteraction(const LiquidTranInteraction &right ) { + *this = right; //use assignment operator to do other work + } + + // Assignment operator + LiquidTranInteraction& LiquidTranInteraction::operator=(const LiquidTranInteraction &right ) + { + if (&right != this) { + m_model = right.m_model; + m_property = right.m_property; + 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; + } + return *this; + } + + + //==================================================================================================================== + LTI_Solvent::LTI_Solvent(TransportPropertyList tp_ind) : + LiquidTranInteraction(tp_ind) + { + m_model = LTI_MODEL_SOLVENT; + } + //==================================================================================================================== + + 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); + + doublereal value = 0.0; + + //if weightings are specified, use those + if (speciesWeight) { + for (int k = 0; k < nsp; k++) { + molefracs[k] = molefracs[k]; + // should be: molefracs[k] = molefracs[k]*speciesWeight[k]; for consistency, but weight(solvent)=1? + } + } + else { + 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++ ){ + //presume that the weighting is set to 1.0 for solvent and 0.0 for everything else. + value += speciesValues[i] * speciesWeight[i]; + if (i == 0) { + AssertTrace(speciesWeight[i] == 1.0); + } else { + AssertTrace(speciesWeight[i] == 0.0); + } + 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_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.0; + + for (int k = 0; k < nsp; k++) { + 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++) { + //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; + } + //==================================================================================================================== + void LTI_Solvent::getMatrixTransProp(DenseMatrix &mat, doublereal* speciesValues) { + mat = (*m_Aij[0]); + } + //==================================================================================================================== + + + doublereal LTI_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_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); + } + } + } + + 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; + } + + + doublereal LTI_Log_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; + + //if weightings are specified, use those + + for (int k = 0; k < nsp; k++) { + molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight(); + } + + 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; + } + + + + + + void LTI_Pairwise_Interaction::setParameters(LiquidTransportParams& trParam) { + int nsp = m_thermo->nSpecies(); + m_diagonals.resize(nsp, 0); + + for (int k = 0; k < nsp; k++) { + Cantera::LiquidTransportData <d = trParam.LTData[k]; + if (ltd.speciesDiffusivity) + m_diagonals[k] = ltd.speciesDiffusivity; + } + } + + doublereal LTI_Pairwise_Interaction::getMixTransProp(doublereal *speciesValues, doublereal *speciesWeight) { + + int nsp = m_thermo->nSpecies(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense." ); + + return value; + } + + + doublereal LTI_Pairwise_Interaction::getMixTransProp(std::vector LTPptrs) { + + int nsp = m_thermo->nSpecies(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense." ); + + return value; + } + + void LTI_Pairwise_Interaction::getMatrixTransProp(DenseMatrix &mat, doublereal *speciesValues ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs ); + + mat.resize(nsp, nsp, 0.0 ); + for (int i = 0; i < nsp; i++ ) + for (int j = 0; j < i; j++ ) + mat(i,j) = mat(j,i) = exp(m_Eij(i,j) / temp ) / m_Dij(i,j); + + for (int i = 0; i < nsp; i++ ) + if (mat(i,i) == 0.0 && m_diagonals[i] ) + mat(i,i) = 1.0 / m_diagonals[i]->getSpeciesTransProp() ; + } + + + void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam ) { + int nsp = m_thermo->nSpecies(); + int nsp2 = nsp*nsp; + //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 nsp2 = nsp*nsp; + 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); + std::vector neutMolIndex(3); + ions_thermo->getDissociationCoeffs(viS,charges,neutMolIndex); + + 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"); + } + + m_ionCondMix = m_ionCondMixModel->getMixTransProp(m_ionCondSpecies); + + MargulesVPSSTP * marg_thermo = dynamic_cast (ions_thermo->neutralMoleculePhase_); + doublereal vol = m_thermo->molarVolume(); + + k = 0; + for (j = 0; j < nsp; j++) { + for (i = 0; i < nsp; i++) { + if (m_mobRatMixModel[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); + } + } + k++; + } + } + + + for (k = 0; k < nsp; k++) { + m_selfDiffMix[k] = m_selfDiffMixModel[k]->getMixTransProp(m_selfDiffSpecies[k] ); + } + + 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); + marg_thermo->getdlnActCoeffdlnN(&dlnActCoeffdlnN[0]); + + xA = neut_molefracs[neutMolIndex[cation[0]]]; + xB = neut_molefracs[neutMolIndex[cation[1]]]; + eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0])); + inv_vP_vM_MutualDiff = (xA*(1+dlnActCoeffdlnN[neutMolIndex[cation[1]]])/m_selfDiffMix[cation[1]]+xB*(1+dlnActCoeffdlnN[neutMolIndex[cation[0]]])/m_selfDiffMix[cation[0]]); + + 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/m_ionCondMix/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/m_ionCondMix/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/m_ionCondMix/vol; + + } + + + doublereal LTI_StokesEinstein::getMixTransProp(doublereal *speciesValues, doublereal *speciesWeight ) { + + int nsp = m_thermo->nSpecies(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense." ); + + return value; + } + + + doublereal LTI_StokesEinstein::getMixTransProp(std::vector LTPptrs ) { + + int nsp = m_thermo->nSpecies(); + doublereal molefracs[nsp]; + m_thermo->getMoleFractions(molefracs); + + doublereal value = 0; + + throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense." ); + + return value; + } + + + + void LTI_StokesEinstein::setParameters(LiquidTransportParams& trParam ) { + int nsp = m_thermo->nSpecies(); + m_viscosity.resize(nsp, 0 ); + m_hydroRadius.resize(nsp, 0 ); + for (int k = 0; k < nsp; k++) { + Cantera::LiquidTransportData <d = trParam.LTData[k]; + m_viscosity[k] = ltd.viscosity; + m_hydroRadius[k] = ltd.hydroRadius; + } + } + + void LTI_StokesEinstein::getMatrixTransProp(DenseMatrix &mat, doublereal *speciesValues ) { + + int nsp = m_thermo->nSpecies(); + doublereal temp = m_thermo->temperature(); + + double *viscSpec = new double(nsp); + double *radiusSpec = new double(nsp); + + for (int k = 0; k < nsp; k++) { + viscSpec[k] = m_viscosity[k]->getSpeciesTransProp() ; + radiusSpec[k] = m_hydroRadius[k]->getSpeciesTransProp() ; + } + + mat.resize(nsp,nsp, 0.0); + for (int i = 0; i < nsp; i++) + for (int j = 0; j < nsp; j++) { + mat(i,j) = (6.0 * Pi * radiusSpec[i] * viscSpec[j] ) / GasConstant / temp; + } + delete radiusSpec; + 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; + } + + + +} //namespace Cantera diff --git a/Cantera/src/transport/LiquidTranInteraction.h b/Cantera/src/transport/LiquidTranInteraction.h new file mode 100644 index 000000000..f7d434c0d --- /dev/null +++ b/Cantera/src/transport/LiquidTranInteraction.h @@ -0,0 +1,685 @@ +/** + * @file LiquidTranInteraction.h + * Header file defining the class LiquidTranInteraction and classes which + * derive from LiquidTranInteraction. + */ +/* + * $Author: christopher.lueth@gmail.com $ + * $Date: 2010-03-30 12:30:39 -0600 (Tue, 30 Mar 2010) $ + * $Revision: 427 $ + * + */ +#ifndef CT_LIQUIDTRANINTERACTION_H +#define CT_LIQUIDTRANINTERACTION_H + +#include "ct_defs.h" +#include "TransportBase.h" +#include "TransportParams.h" +#include "LiquidTransportData.h" +#include "xml.h" +#include "XML_Writer.h" + +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 + * -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, + LTI_MODEL_SOLVENT, + LTI_MODEL_MOLEFRACS, + LTI_MODEL_MASSFRACS, + LTI_MODEL_LOG_MOLEFRACS, + LTI_MODEL_PAIRWISE_INTERACTION, + LTI_MODEL_STEFANMAXWELL_PPN, + LTI_MODEL_STOKES_EINSTEIN, + LTI_MODEL_MOLEFRACS_EXPT + }; + + + //! Base class to handle transport property evaluation in a mixture. + /*! + * In a mixture, the mixture transport properties will generally depend on + * the contributions of each of the pure species transport properties. + * Many composition dependencies are possible. This class, + * LiquidTranInteraction, is designed to be a base class for the + * implementation of various models for the mixing of pure species + * transport properties. + * + * There are two very broad types of transport properties to consider. + * First, there are properties for which a mixture value can be + * obtained through some mixing rule. These are obtained using the + * method getMixTransProp(). Viscosity is typical of this. + * Second there are properties for which a matrix of properties may + * @param tp_ind + * exist. This matrix of properties is obtained from the method + * getMatrixTransProp(). Diffusion coefficients are of this type. + * Subclasses should implement the appropriate one or both of + * these methods. + * + */ + class LiquidTranInteraction { + + public: + //! Constructor + /** + * @param tp_ind Index indicating transport property type (i.e. viscosity) + */ + LiquidTranInteraction(TransportPropertyList tp_ind = TP_UNKNOWN); + + //! Copy constructor + LiquidTranInteraction(const LiquidTranInteraction &right); + + //! Assignment operator + LiquidTranInteraction& operator=(const LiquidTranInteraction &right ); + + //! destructor + virtual ~LiquidTranInteraction(); + + //! initialize LiquidTranInteraction objects with thermo and XML node + /** + * @param compModelNode \verbatim \endverbatim XML node + * @param thermo Pointer to thermo object + */ + virtual void init(const XML_Node &compModelNode = 0, + thermo_t* thermo = 0 ); + + virtual void setParameters( LiquidTransportParams& trParam ) { ; } + + //! Return the mixture transport property value. + //! (Must be implemented in subclasses.) + virtual doublereal getMixTransProp(doublereal* speciesValues, doublereal *weightSpecies = 0 ) { + throw NotImplemented("LiquidTranInteraction::getMixTransProp"); + } + + virtual doublereal getMixTransProp(std::vector LTPptrs ) { + throw NotImplemented("LiquidTranInteraction::getMixTransProp"); + } + + virtual void getMatrixTransProp(DenseMatrix &mat, doublereal* speciesValues = 0 ) { + //mat = m_Dij; + throw NotImplemented("LiquidTranInteraction::getMixTransProp"); + } + + protected: + //! Model for species interaction effects + //! Takes enum LiquidTranMixingModel + LiquidTranMixingModel m_model; + + //! enum indicating what property this is (i.e viscosity) + TransportPropertyList m_property; + + //! pointer to thermo object to get current temperature + thermo_t* m_thermo; + + //LiquidTransportParams* m_trParam; + + //! 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 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; + + //! Matrix of interactions + DenseMatrix m_Dij; + }; + + class LTI_Solvent : public LiquidTranInteraction { + + public: + LTI_Solvent(TransportPropertyList tp_ind = TP_UNKNOWN); + + //! Copy constructor + // LTI_Solvent( const LTI_Solvent &right ); + + //! Assignment operator + // LTI_Solvent& operator=( const LTI_Solvent &right ); + + virtual ~LTI_Solvent( ) { } + + //! 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) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them: Not implemented for this mixing rule. + */ + void getMatrixTransProp(DenseMatrix &mat, doublereal* speciesValues = 0 ); + + protected: + + }; + + //! 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 : public LiquidTranInteraction { + + public: + LTI_MoleFracs( TransportPropertyList tp_ind = TP_UNKNOWN ) : + LiquidTranInteraction( tp_ind ) + { + m_model = LTI_MODEL_MOLEFRACS; + } + + + //! Copy constructor + // LTI_MoleFracs( const LTI_MoleFracs &right ); + + //! Assignment operator + // LTI_MoleFracs& operator=( const LTI_MoleFracs &right ); + + virtual ~LTI_MoleFracs( ) { } + + //! 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 ) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them: Not Implemented for this Mixing rule; + */ + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } + + protected: + + }; + + + //! Simple mass fraction weighting of transport properties + /*! + * This model weights the transport property by the mass + * fractions. + * The overall formula for the mixture viscosity is + * + * \f[ + * \eta_{mix} = \sum_i Y_i \eta_i + * + \sum_i \sum_j Y_i Y_j A_{i,j} + * \f]. + */ + class LTI_MassFracs : public LiquidTranInteraction { + + public: + + LTI_MassFracs( TransportPropertyList tp_ind = TP_UNKNOWN ) : + LiquidTranInteraction( tp_ind ) + { + m_model = LTI_MODEL_MASSFRACS; + } + + + //! Copy constructor + // LTI_MassFracs( const LTI_MassFracs &right ); + + //! Assignment operator + // LTI_MassFracs& operator=( const LTI_MassFracs &right ); + + virtual ~LTI_MassFracs( ) { } + + //! 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 ) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them: Not implemented for this mixing rule. + */ + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } + + protected: + + }; + + + //! Mixing rule using logarithms of the mole fractions + /** + * This model is based on the idea that liquid molecules are + * generally interacting with some energy and entropy of interaction. + * For transport properties that depend on these energies of + * interaction, the mixture transport property can be written + * in terms of its logarithm + * + * \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]. + * + * These additional interaction terms multiply the mixture property by + * \f[ \exp( \sum_{i} \sum_{j} X_i X_j ( S_{i,j} + E_{i,j} / T ) ) \f] + * so that the self-interaction terms \f$ S_{i,j} \f$ and + * \f$ E_{i,j} \f$ should be zero. + * + * Note that the energies and entropies of interaction should be + * a function of the composition themselves, but this is not yet + * implemented. (We might follow the input of Margules model + * thermodynamic data for the purpose of implementing this.) + * + * Sample input for this method is + * \verbatim + * + * + * + * + * + * -1.0e3 + * 80.0e-5 + * + * + * + * + * \endverbatim + */ + class LTI_Log_MoleFracs : public LiquidTranInteraction { + + public: + LTI_Log_MoleFracs( TransportPropertyList tp_ind = TP_UNKNOWN ) : + LiquidTranInteraction( tp_ind ) + { + m_model = LTI_MODEL_LOG_MOLEFRACS; + } + + + //! Copy constructor + // LTI_Log_MoleFracs( const LTI_Log_MoleFracs &right ); + + //! Assignment operator + // LTI_Log_MoleFracs& operator=( const LTI_Log_MoleFracs &right ); + + virtual ~LTI_Log_MoleFracs( ) { } + + //! 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 ) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them: Not implemented for this mixing rule. + */ + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = m_Eij; } + + protected: + + }; + + + //! Transport properties that act like pairwise interactions + //! as in binary diffusion coefficients. + /** + * This class holds parameters for transport properties expressed + * as a matrix of pairwise interaction parameters. + * Input can be provided for constant or Arrhenius forms of the + * separate parameters. + * + * Sample input for this method is + * \verbatim + * + * + * + * + * 1.0e-8 + * 24.0e6 + * + * + * + * + * \endverbatim + * + */ + class LTI_Pairwise_Interaction : public LiquidTranInteraction { + + public: + LTI_Pairwise_Interaction( TransportPropertyList tp_ind = TP_UNKNOWN ) : + LiquidTranInteraction( tp_ind ) + { + m_model = LTI_MODEL_PAIRWISE_INTERACTION; + } + + + //! Copy constructor + // LTI_Pairwise_Interaction( const LTI_Pairwise_Interaction &right ); + + //! Assignment operator + // LTI_Pairwise_Interaction& operator=( const LTI_Pairwise_Interaction &right ); + + virtual ~LTI_Pairwise_Interaction( ) { } + + 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 ) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them + */ + 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 class is only valid for a common anion mixture of two + //! salts with cations of equal charge. Hence the name _PPN. + /** + * + * This class requres you specify + * + * 1 - ion conductivity + * + * 2 - mobility ratio of the two cations (set all other ratios to zero) + * + * 3 - Self diffusion coefficients of the cations (set others to zero) + * is used to calculate the "mutual diffusion coefficient". The + * approximation needed to do so requires the cations have equal charge. + * + * We than calculate the Stefan Maxwell Diffusion Coefficients by + * \f[ + * \frac{1}{D_{12}} = (1-\epsilon X_A)(1+\epsilon X_B) + * \frac{\nu_- + \nu_+}{\nu_-\nu_+^2D} + * + \frac{z_-z_+ F^2}{\kappa V R T} + * \f] + * \f[ + * \frac{1}{D_{12}} = -\epsilon X_B(1-\epsilon X_A) + * \frac{\nu_- + \nu_+}{\nu_-^2\nu_+D} + * - \frac{z_-z_+ F^2}{\kappa V R T} + * \f] + * \f[ + * \frac{1}{D_{23}} = \epsilon X_A(1+\epsilon X_B) + * \frac{\nu_- + \nu_+}{\nu_-^2\nu_+D} + * - \frac{z_-z_+ F^2}{\kappa V R T} + * \f] + * where F is Faraday's constant, RT is the gas constant times the + * tempurature, and V is the molar volume (basis is moles of ions) that is + * calculated by the thermophase member. X_A and X_B are the mole fractions + * of the salts composed of cation(1) and cation(2), respectively, that share + * a common anion(3). \f$\nu_{+,-}\f$ are the stoichiometric coefficients in + * the dissociation reaction of the salts to the ions with charges of + * \f$z_{+,-}\f$. Assuming that the cations have equal charge, the "mutual + * diffusion coefficient" is calculated using the cation self diffusion + * coefficients. + * \f[ + * \frac{1}{\nu_-\nu_+D} = \left(1+\frac{\partial \gamma_B}{\partial N_B} + * \right)\frac{X_A}{D_2^*}+\left(1+\frac{\partial \gamma_A}{\partial N_A} + * \right)\frac{X_B}{D_1^*} + * \f] + * where the self diffusion coefficients, \f$D_i^*\f$, are temperature and + * composition parameterized inputs and the derivative of the activity + * coefficient, \f$\frac{\partial \gamma_B}{\partial N_B}\f$, is calculated + * by the thermophase member using the excess enthalpy and entropy upon mixing. + * + * Finally, the deviation of the transferrence numbers from ideality, + * \f$\epsilon\f$, is calculated from the mobility ratio of the cations. + * \f[ + * \epsilon = \frac{1-b_2/b_1}{X_A+X_Bb_2/b_1} + * \f] + * Where \f$b_i\f$ are the mobilities of the two cations. Everywhere, + * cation 1 corresponds with salt A and cation 2 with salt B. + * + * 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 ) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them + */ + void getMatrixTransProp( DenseMatrix &mat, 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_selfDiffMixModel; + vector_fp m_selfDiffMix; + std::vector m_selfDiffSpecies; + }; + + + class LTI_StokesEinstein : public LiquidTranInteraction { + + public: + LTI_StokesEinstein( TransportPropertyList tp_ind = TP_UNKNOWN ) : + LiquidTranInteraction( tp_ind ) + { + m_model = LTI_MODEL_STOKES_EINSTEIN; + } + + + //! Copy constructor + // LTI_StokesEinstein( const LTI_StokesEinstein &right ); + + //! Assignment operator + // LTI_StokesEinstein& operator=( const LTI_StokesEinstein &right ); + + virtual ~LTI_StokesEinstein( ) { } + + 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 ) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them + */ + void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) ; + protected: + + std::vector m_viscosity; + std::vector m_hydroRadius; + + }; + + //! 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 ) ; + + //! Return the matrix of binary interaction parameters. + /** + * Takes the proper mixing rule for the binary interaction parameters + * and calculates them: Not Implemented for this mixing rule + */ + 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: + + }; + +} + +#endif diff --git a/Cantera/src/transport/LiquidTransportParams.cpp b/Cantera/src/transport/LiquidTransportParams.cpp index f0c587060..73323b97e 100644 --- a/Cantera/src/transport/LiquidTransportParams.cpp +++ b/Cantera/src/transport/LiquidTransportParams.cpp @@ -19,213 +19,30 @@ using namespace std; namespace Cantera { - /** - * Exception thrown if an error is encountered while reading the - * transport database. - */ - class LTPError : public CanteraError { - public: - LTPError( std::string msg ) - : CanteraError("LTPspecies", - "error parsing transport data: " - + msg + "\n") {} - }; - - /** - * Exception thrown if an error is encountered while reading the - * transport database. - */ + //==================================================================================================================== + //! Exception thrown if an error is encountered while reading the transport database. class LTPmodelError : public CanteraError { public: - LTPmodelError( std::string msg ) - : CanteraError("LTPspecies", - "error parsing transport data: " - + msg + "\n") {} + LTPmodelError(std::string msg) : + CanteraError("LTPspecies", "error parsing transport data: " + msg + "\n") + { + } }; - - - // Constructor - /* - * @param tp_ind Index indicating transport property type (i.e. viscosity) - */ - LiquidTranInteraction::LiquidTranInteraction( TransportPropertyList tp_ind ) : - m_model(LTI_MODEL_NOTSET), - m_property(tp_ind) - { - } - - 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_Dij.resize( nsp, nsp, 0.0 ); - m_Eij.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; - - int num = compModelNode.nChildren(); - for (int iChild = 0; iChild < num; iChild++) { - XML_Node &xmlChild = compModelNode.child(iChild); - std::string nodeName = lowercase( xmlChild.name() ); - if ( nodeName != "interaction" ) - throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", - "expected element and got <" + nodeName + ">" ); - speciesA = xmlChild.attrib("speciesA"); - speciesB = xmlChild.attrib("speciesB"); - int iSpecies = m_thermo->speciesIndex( speciesA ); - if ( iSpecies < 0 ) - throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", - "Unknown species " + speciesA ); - int jSpecies = m_thermo->speciesIndex( speciesB ); - if ( jSpecies < 0 ) - throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", - "Unknown species " + speciesB ); - /* 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" ); - m_Eij(iSpecies,jSpecies) /= GasConstant; - m_Eij(jSpecies,iSpecies) = m_Eij(iSpecies,jSpecies) ; - } - - 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) ; - }*/ - - if ( xmlChild.hasChild( "Dij" ) ) { - m_Dij(iSpecies,jSpecies) = getFloat( xmlChild, "Dij", "toSI" ); - m_Dij(jSpecies,iSpecies) = m_Dij(iSpecies,jSpecies) ; - } - } - } - - // Copy constructor - LiquidTranInteraction::LiquidTranInteraction( const LiquidTranInteraction &right ) { - *this = right; //use assignment operator to do other work - } - - // Assignment operator - LiquidTranInteraction& LiquidTranInteraction::operator=( const LiquidTranInteraction &right ) - { - if (&right != this) { - m_model = right.m_model; - m_property = right.m_property; - 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; - } - return *this; - } - //==================================================================================================================== LiquidTransportParams::LiquidTransportParams() : - viscosity(0), ionConductivity(0), thermalCond(0), speciesDiffusivity(0), electCond(0), hydroRadius(0), model_viscosity(LTI_MODEL_NOTSET), - model_speciesDiffusivity(LTI_MODEL_NOTSET), model_hydroradius(LTI_MODEL_NOTSET) + TransportParams(), + LTData(0), + viscosity(0), + ionConductivity(0), + mobilityRatio(0), + selfDiffusion(0), + thermalCond(0), + speciesDiffusivity(0), + electCond(0), + hydroRadius(0), + model_viscosity(LTI_MODEL_NOTSET), + model_speciesDiffusivity(LTI_MODEL_NOTSET), + model_hydroradius(LTI_MODEL_NOTSET) { } @@ -239,15 +56,21 @@ namespace Cantera { delete electCond; delete hydroRadius; } - //==================================================================================================================== LiquidTransportParams::LiquidTransportParams(const LiquidTransportParams &right) : - viscosity(0), thermalCond(0), speciesDiffusivity(0), electCond(0), hydroRadius(0), model_viscosity(LTI_MODEL_NOTSET), - model_speciesDiffusivity(LTI_MODEL_NOTSET), model_hydroradius(LTI_MODEL_NOTSET) + TransportParams(), + LTData(0), + viscosity(0), + thermalCond(0), + speciesDiffusivity(0), + electCond(0), + hydroRadius(0), + model_viscosity(LTI_MODEL_NOTSET), + model_speciesDiffusivity(LTI_MODEL_NOTSET), + model_hydroradius(LTI_MODEL_NOTSET) { - throw CanteraError("LiquidTransportParams(const LiquidTransportParams &right)","not implemented"); + operator=(right); } - //==================================================================================================================== LiquidTransportParams& LiquidTransportParams::operator=(const LiquidTransportParams & right) @@ -256,631 +79,50 @@ namespace Cantera { return *this; } - throw CanteraError("LiquidTransportParams(const LiquidTransportParams &right)","not implemented"); + LTData = right.LTData; + + delete viscosity; + if (right.viscosity) { + viscosity = new LiquidTranInteraction(*(right.viscosity)); + } + delete ionConductivity; + if (right.ionConductivity) { + ionConductivity = new LiquidTranInteraction(*(right.ionConductivity)); + } + deepStdVectorPointerCopy(right.mobilityRatio, mobilityRatio); + deepStdVectorPointerCopy(right.selfDiffusion, selfDiffusion); + + delete thermalCond; + if (right.thermalCond) { + thermalCond = new LiquidTranInteraction(*(right.thermalCond)); + } + delete speciesDiffusivity; + if (right.speciesDiffusivity) { + speciesDiffusivity = new LiquidTranInteraction(*(right.speciesDiffusivity)); + } + + delete electCond; + if (right.electCond) { + electCond = new LiquidTranInteraction(*(right.electCond)); + } + delete hydroRadius; + if (right.hydroRadius) { + hydroRadius = new LiquidTranInteraction(*(right.hydroRadius)); + } + model_viscosity = right.model_viscosity; + model_ionConductivity = right.model_ionConductivity; + deepStdVectorPointerCopy(right.model_mobilityRatio, model_mobilityRatio); + deepStdVectorPointerCopy(right.model_selfDiffusion, model_selfDiffusion); + thermalCond_Aij = right.thermalCond_Aij; + model_speciesDiffusivity = right.model_speciesDiffusivity; + diff_Dij = right.diff_Dij; + model_hydroradius = right.model_hydroradius; + radius_Aij = right.radius_Aij; + + throw CanteraError("LiquidTransportParams(const LiquidTransportParams &right)", "not tested"); + return *this; - } - //==================================================================================================================== - - - 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); - - doublereal value = 0; - - //if weightings are specified, use those - if ( speciesWeight ) { - for ( int k = 0; k < nsp; k++) { - molefracs[k] = molefracs[k]; - // should be: molefracs[k] = molefracs[k]*speciesWeight[k]; for consistency, but weight(solvent)=1? - } - } - else { - 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++ ){ - //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; - } - - - 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++) { - 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++ ){ - //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; - } - - - - - - doublereal LTI_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_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); - } - } - } - - 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; - } - - - doublereal LTI_Log_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; - - //if weightings are specified, use those - - for ( int k = 0; k < nsp; k++) { - molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight( ); - } - - 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; - } - - - - - - void LTI_Pairwise_Interaction::setParameters(LiquidTransportParams& trParam) { - int nsp = m_thermo->nSpecies(); - m_diagonals.resize(nsp, 0); - - for (int k = 0; k < nsp; k++) { - Cantera::LiquidTransportData <d = trParam.LTData[k]; - if (ltd.speciesDiffusivity) - m_diagonals[k] = ltd.speciesDiffusivity; - } - } - - doublereal LTI_Pairwise_Interaction::getMixTransProp(doublereal *speciesValues, doublereal *speciesWeight) { - - int nsp = m_thermo->nSpecies(); - doublereal molefracs[nsp]; - m_thermo->getMoleFractions(molefracs); - - doublereal value = 0; - - throw LTPmodelError( "Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense." ); - - return value; - } - - - doublereal LTI_Pairwise_Interaction::getMixTransProp(std::vector LTPptrs) { - - int nsp = m_thermo->nSpecies(); - doublereal molefracs[nsp]; - m_thermo->getMoleFractions(molefracs); - - doublereal value = 0; - - throw LTPmodelError( "Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense." ); - - return value; - } - - void LTI_Pairwise_Interaction::getMatrixTransProp( DenseMatrix &mat, doublereal *speciesValues ) { - - int nsp = m_thermo->nSpecies(); - doublereal temp = m_thermo->temperature(); - doublereal molefracs[nsp]; - m_thermo->getMoleFractions( molefracs ); - - mat.resize( nsp, nsp, 0.0 ); - for ( int i = 0; i < nsp; i++ ) - for ( int j = 0; j < i; j++ ) - mat(i,j) = mat(j,i) = exp( m_Eij(i,j) / temp ) / m_Dij(i,j); - - for ( int i = 0; i < nsp; i++ ) - if ( mat(i,i) == 0.0 && m_diagonals[i] ) - mat(i,i) = 1.0 / m_diagonals[i]->getSpeciesTransProp() ; - } - - - void LTI_StefanMaxwell_PPN::setParameters( LiquidTransportParams& trParam ) { - int nsp = m_thermo->nSpecies(); - int nsp2 = nsp*nsp; - //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 nsp2 = nsp*nsp; - 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); - std::vector neutMolIndex(3); - ions_thermo->getDissociationCoeffs(viS,charges,neutMolIndex); - - 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"); - } - - m_ionCondMix = m_ionCondMixModel->getMixTransProp(m_ionCondSpecies); - - MargulesVPSSTP * marg_thermo = dynamic_cast (ions_thermo->neutralMoleculePhase_); - doublereal vol = m_thermo->molarVolume(); - - k = 0; - for (j = 0; j < nsp; j++) { - for (i = 0; i < nsp; i++) { - if (m_mobRatMixModel[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); - } - } - k++; - } - } - - - for (k = 0; k < nsp; k++) { - m_selfDiffMix[k] = m_selfDiffMixModel[k]->getMixTransProp( m_selfDiffSpecies[k] ); - } - - 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); - marg_thermo->getdlnActCoeffdlnN(&dlnActCoeffdlnN[0]); - - xA = neut_molefracs[neutMolIndex[cation[0]]]; - xB = neut_molefracs[neutMolIndex[cation[1]]]; - eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0])); - inv_vP_vM_MutualDiff = (xA*(1+dlnActCoeffdlnN[neutMolIndex[cation[1]]])/m_selfDiffMix[cation[1]]+xB*(1+dlnActCoeffdlnN[neutMolIndex[cation[0]]])/m_selfDiffMix[cation[0]]); - - 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/m_ionCondMix/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/m_ionCondMix/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/m_ionCondMix/vol; - - } - - - doublereal LTI_StokesEinstein::getMixTransProp( doublereal *speciesValues, doublereal *speciesWeight ) { - - int nsp = m_thermo->nSpecies(); - doublereal molefracs[nsp]; - m_thermo->getMoleFractions(molefracs); - - doublereal value = 0; - - throw LTPmodelError( "Calling LTI_StokesEinstein::getMixTransProp does not make sense." ); - - return value; - } - - - doublereal LTI_StokesEinstein::getMixTransProp( std::vector LTPptrs ) { - - int nsp = m_thermo->nSpecies(); - doublereal molefracs[nsp]; - m_thermo->getMoleFractions(molefracs); - - doublereal value = 0; - - throw LTPmodelError( "Calling LTI_StokesEinstein::getMixTransProp does not make sense." ); - - return value; - } - - - - void LTI_StokesEinstein::setParameters( LiquidTransportParams& trParam ) { - int nsp = m_thermo->nSpecies(); - m_viscosity.resize( nsp, 0 ); - m_hydroRadius.resize( nsp, 0 ); - for (int k = 0; k < nsp; k++) { - Cantera::LiquidTransportData <d = trParam.LTData[k]; - m_viscosity[k] = ltd.viscosity; - m_hydroRadius[k] = ltd.hydroRadius; - } - } - - void LTI_StokesEinstein::getMatrixTransProp( DenseMatrix &mat, doublereal *speciesValues ) { - - int nsp = m_thermo->nSpecies(); - doublereal temp = m_thermo->temperature(); - - double *viscSpec = new double(nsp); - double *radiusSpec = new double(nsp); - - for (int k = 0; k < nsp; k++) { - viscSpec[k] = m_viscosity[k]->getSpeciesTransProp() ; - radiusSpec[k] = m_hydroRadius[k]->getSpeciesTransProp() ; - } - - mat.resize(nsp,nsp, 0.0); - for (int i = 0; i < nsp; i++) - for (int j = 0; j < nsp; j++) { - mat(i,j) = ( 6.0 * Pi * radiusSpec[i] * viscSpec[j] ) / GasConstant / temp; - } - delete radiusSpec; - 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; - } - - - } //namespace Cantera diff --git a/Cantera/src/transport/LiquidTransportParams.h b/Cantera/src/transport/LiquidTransportParams.h index ba68fcd2a..eaf1a795d 100644 --- a/Cantera/src/transport/LiquidTransportParams.h +++ b/Cantera/src/transport/LiquidTransportParams.h @@ -13,25 +13,21 @@ #ifndef CT_LIQUIDTRANSPORTPARAMS_H #define CT_LIQUIDTRANSPORTPARAMS_H -#include - #include "ct_defs.h" #include "TransportBase.h" #include "TransportParams.h" #include "LiquidTransportData.h" +#include "LiquidTranInteraction.h" #include "xml.h" #include "XML_Writer.h" + +#include + + namespace Cantera { - class NotImplemented : public CanteraError { - public: - NotImplemented(std::string method) : CanteraError("Transport", - "\n\n**** Method "+method+" not implemented. ****\n" - "(Did you forget to specify a transport model?)\n\n") {} - }; - //! Composition dependence type for liquid mixture transport properties /*! @@ -89,137 +85,45 @@ namespace Cantera { * */ - enum LiquidTranMixingModel { - LTI_MODEL_NOTSET=-1, - LTI_MODEL_NONE, - LTI_MODEL_SOLVENT, - LTI_MODEL_MOLEFRACS, - LTI_MODEL_MASSFRACS, - LTI_MODEL_LOG_MOLEFRACS, - LTI_MODEL_PAIRWISE_INTERACTION, - LTI_MODEL_STEFANMAXWELL_PPN, - LTI_MODEL_STOKES_EINSTEIN, - LTI_MODEL_MOLEFRACS_EXPT - }; - - //! Base class to handle transport property evaluation in a mixture. - /** - * In a mixture, the mixture transport properties will generally depend on - * the contributions of each of the pure species transport properties. - * Many composition dependencies are possible. This class, - * LiquidTranInteraction, is designed to be a base class for the - * implementation of various models for the mixing of pure species - * transport properties. - * - * There are two very broad types of transport properties to consider. - * First, there are properties for which a mixture value can be - * obtained through some mixing rule. These are obtained using the - * method getMixTransProp(). Viscosity is typical of this. - * Second there are properties for which a matrix of properties may - * @param tp_ind - * exist. This matrix of properties is obtained from the method - * getMatrixTransProp(). Diffusion coefficients are of this type. - * Subclasses should implement the appropriate one or both of - * these methods. - * - */ - class LiquidTranInteraction { - - public: - //! Constructor - /** - * @param tp_ind Index indicating transport property type (i.e. viscosity) - */ - LiquidTranInteraction( TransportPropertyList tp_ind = TP_UNKNOWN ); - - //! Copy constructor - LiquidTranInteraction( const LiquidTranInteraction &right ); - - //! Assignment operator - LiquidTranInteraction& operator=( const LiquidTranInteraction &right ); - - //! destructor - virtual ~LiquidTranInteraction(); - //! initialize LiquidTranInteraction objects with thermo and XML node - /** - * @param compModelNode \verbatim \endverbatim XML node - * @param thermo Pointer to thermo object - */ - virtual void init( const XML_Node &compModelNode = 0, - thermo_t* thermo = 0 ); - - virtual void setParameters( LiquidTransportParams& trParam ) { ; } - - //! Return the mixture transport property value. - //! (Must be implemented in subclasses.) - virtual doublereal getMixTransProp( doublereal* speciesValues, doublereal *weightSpecies = 0 ) { - throw NotImplemented("LiquidTranInteraction::getMixTransProp"); - } - - virtual doublereal getMixTransProp( std::vector LTPptrs ) { - throw NotImplemented("LiquidTranInteraction::getMixTransProp"); - } - - virtual void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { - //mat = m_Dij; - throw NotImplemented("LiquidTranInteraction::getMixTransProp"); - } - - protected: - //! Model for species interaction effects - //! Takes enum LiquidTranMixingModel - LiquidTranMixingModel m_model; - - //! enum indicating what property this is (i.e viscosity) - TransportPropertyList m_property; - - //! pointer to thermo object to get current temperature - thermo_t* m_thermo; - - //LiquidTransportParams* m_trParam; - - //! 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 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; - - }; - - /** - * Class LiquidTransportParams holds transport model parameters - * relevant to transport in mixtures. - * Used by TransportFactory. + + //! Class LiquidTransportParams holds transport model parameters + //! relevant to transport in mixtures. + /*! + * This class is used by TransportFactory to initialize transport objects. */ class LiquidTransportParams : public TransportParams { public: + //! Constructor LiquidTransportParams(); + + //! Destructor ~LiquidTransportParams(); + + //! Copy constructor + /*! + * @param right Object to be copied + */ LiquidTransportParams(const LiquidTransportParams &right); + + //! Assignment operator + /*! + * @param right Object to be copied + */ LiquidTransportParams & operator=(const LiquidTransportParams &right); //! Species transport parameters std::vector LTData; + //! Object that specifies the viscosity interaction for the mixture LiquidTranInteraction* viscosity; + + //! Object that specifes the ionic Conductivity of the mixture LiquidTranInteraction* ionConductivity; + std::vector mobilityRatio; std::vector selfDiffusion; LiquidTranInteraction* thermalCond; @@ -279,501 +183,8 @@ namespace Cantera { DenseMatrix radius_Aij; }; - - class LTI_Solvent; - class LTI_MoleFracs; - 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: - LTI_Solvent( TransportPropertyList tp_ind = TP_UNKNOWN ) : - LiquidTranInteraction( tp_ind ) - { - m_model = LTI_MODEL_SOLVENT; - } - - //! Copy constructor - // LTI_Solvent( const LTI_Solvent &right ); - - //! Assignment operator - // LTI_Solvent& operator=( const LTI_Solvent &right ); - - virtual ~LTI_Solvent( ) { } - - //! 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them: Not implemented for this mixing rule. - */ - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } - - protected: - - }; - - //! 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 : public LiquidTranInteraction { - - public: - LTI_MoleFracs( TransportPropertyList tp_ind = TP_UNKNOWN ) : - LiquidTranInteraction( tp_ind ) - { - m_model = LTI_MODEL_MOLEFRACS; - } - - - //! Copy constructor - // LTI_MoleFracs( const LTI_MoleFracs &right ); - - //! Assignment operator - // LTI_MoleFracs& operator=( const LTI_MoleFracs &right ); - - virtual ~LTI_MoleFracs( ) { } - - //! 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them: Not Implemented for this Mixing rule; - */ - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } - - protected: - - }; - - - //! Simple mass fraction weighting of transport properties - /** - * This model weights the transport property by the mass - * fractions. - * The overall formula for the mixture viscosity is - * - * \f[ \eta_{mix} = \sum_i Y_i \eta_i - * + \sum_i \sum_j Y_i Y_j A_{i,j} \f]. - */ - class LTI_MassFracs : public LiquidTranInteraction { - - public: - LTI_MassFracs( TransportPropertyList tp_ind = TP_UNKNOWN ) : - LiquidTranInteraction( tp_ind ) - { - m_model = LTI_MODEL_MASSFRACS; - } - - - //! Copy constructor - // LTI_MassFracs( const LTI_MassFracs &right ); - - //! Assignment operator - // LTI_MassFracs& operator=( const LTI_MassFracs &right ); - - virtual ~LTI_MassFracs( ) { } - - //! 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them: Not implemented for this mixing rule. - */ - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = (*m_Aij[0]); } - - protected: - - }; - - - //! Mixing rule using logarithms of the mole fractions - /** - * This model is based on the idea that liquid molecules are - * generally interacting with some energy and entropy of interaction. - * For transport properties that depend on these energies of - * interaction, the mixture transport property can be written - * in terms of its logarithm - * - * \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]. - * - * These additional interaction terms multiply the mixture property by - * \f[ \exp( \sum_{i} \sum_{j} X_i X_j ( S_{i,j} + E_{i,j} / T ) ) \f] - * so that the self-interaction terms \f$ S_{i,j} \f$ and - * \f$ E_{i,j} \f$ should be zero. - * - * Note that the energies and entropies of interaction should be - * a function of the composition themselves, but this is not yet - * implemented. (We might follow the input of Margules model - * thermodynamic data for the purpose of implementing this.) - * - * Sample input for this method is - * \verbatim - * - * - * - * - * - * -1.0e3 - * 80.0e-5 - * - * - * - * - * \endverbatim - */ - class LTI_Log_MoleFracs : public LiquidTranInteraction { - - public: - LTI_Log_MoleFracs( TransportPropertyList tp_ind = TP_UNKNOWN ) : - LiquidTranInteraction( tp_ind ) - { - m_model = LTI_MODEL_LOG_MOLEFRACS; - } - - - //! Copy constructor - // LTI_Log_MoleFracs( const LTI_Log_MoleFracs &right ); - - //! Assignment operator - // LTI_Log_MoleFracs& operator=( const LTI_Log_MoleFracs &right ); - - virtual ~LTI_Log_MoleFracs( ) { } - - //! 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them: Not implemented for this mixing rule. - */ - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) { mat = m_Eij; } - - protected: - - }; - - - //! Transport properties that act like pairwise interactions - //! as in binary diffusion coefficients. - /** - * This class holds parameters for transport properties expressed - * as a matrix of pairwise interaction parameters. - * Input can be provided for constant or Arrhenius forms of the - * separate parameters. - * - * Sample input for this method is - * \verbatim - * - * - * - * - * 1.0e-8 - * 24.0e6 - * - * - * - * - * \endverbatim - * - */ - class LTI_Pairwise_Interaction : public LiquidTranInteraction { - - public: - LTI_Pairwise_Interaction( TransportPropertyList tp_ind = TP_UNKNOWN ) : - LiquidTranInteraction( tp_ind ) - { - m_model = LTI_MODEL_PAIRWISE_INTERACTION; - } - - - //! Copy constructor - // LTI_Pairwise_Interaction( const LTI_Pairwise_Interaction &right ); - - //! Assignment operator - // LTI_Pairwise_Interaction& operator=( const LTI_Pairwise_Interaction &right ); - - virtual ~LTI_Pairwise_Interaction( ) { } - - 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them - */ - 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 class is only valid for a common anion mixture of two - //! salts with cations of equal charge. Hence the name _PPN. - /** - * - * This class requres you specify - * - * 1 - ion conductivity - * - * 2 - mobility ratio of the two cations (set all other ratios to zero) - * - * 3 - Self diffusion coefficients of the cations (set others to zero) - * is used to calculate the "mutual diffusion coefficient". The - * approximation needed to do so requires the cations have equal charge. - * - * We than calculate the Stefan Maxwell Diffusion Coefficients by - * \f[ - * \frac{1}{D_{12}} = (1-\epsilon X_A)(1+\epsilon X_B) - * \frac{\nu_- + \nu_+}{\nu_-\nu_+^2D} - * + \frac{z_-z_+ F^2}{\kappa V R T} - * \f] - * \f[ - * \frac{1}{D_{12}} = -\epsilon X_B(1-\epsilon X_A) - * \frac{\nu_- + \nu_+}{\nu_-^2\nu_+D} - * - \frac{z_-z_+ F^2}{\kappa V R T} - * \f] - * \f[ - * \frac{1}{D_{23}} = \epsilon X_A(1+\epsilon X_B) - * \frac{\nu_- + \nu_+}{\nu_-^2\nu_+D} - * - \frac{z_-z_+ F^2}{\kappa V R T} - * \f] - * where F is Faraday's constant, RT is the gas constant times the - * tempurature, and V is the molar volume (basis is moles of ions) that is - * calculated by the thermophase member. X_A and X_B are the mole fractions - * of the salts composed of cation(1) and cation(2), respectively, that share - * a common anion(3). \f$\nu_{+,-}\f$ are the stoichiometric coefficients in - * the dissociation reaction of the salts to the ions with charges of - * \f$z_{+,-}\f$. Assuming that the cations have equal charge, the "mutual - * diffusion coefficient" is calculated using the cation self diffusion - * coefficients. - * \f[ - * \frac{1}{\nu_-\nu_+D} = \left(1+\frac{\partial \gamma_B}{\partial N_B} - * \right)\frac{X_A}{D_2^*}+\left(1+\frac{\partial \gamma_A}{\partial N_A} - * \right)\frac{X_B}{D_1^*} - * \f] - * where the self diffusion coefficients, \f$D_i^*\f$, are temperature and - * composition parameterized inputs and the derivative of the activity - * coefficient, \f$\frac{\partial \gamma_B}{\partial N_B}\f$, is calculated - * by the thermophase member using the excess enthalpy and entropy upon mixing. - * - * Finally, the deviation of the transferrence numbers from ideality, - * \f$\epsilon\f$, is calculated from the mobility ratio of the cations. - * \f[ - * \epsilon = \frac{1-b_2/b_1}{X_A+X_Bb_2/b_1} - * \f] - * Where \f$b_i\f$ are the mobilities of the two cations. Everywhere, - * cation 1 corresponds with salt A and cation 2 with salt B. - * - * 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them - */ - void getMatrixTransProp( DenseMatrix &mat, 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_selfDiffMixModel; - vector_fp m_selfDiffMix; - std::vector m_selfDiffSpecies; - }; - - - class LTI_StokesEinstein : public LiquidTranInteraction { - - public: - LTI_StokesEinstein( TransportPropertyList tp_ind = TP_UNKNOWN ) : - LiquidTranInteraction( tp_ind ) - { - m_model = LTI_MODEL_STOKES_EINSTEIN; - } - - - //! Copy constructor - // LTI_StokesEinstein( const LTI_StokesEinstein &right ); - - //! Assignment operator - // LTI_StokesEinstein& operator=( const LTI_StokesEinstein &right ); - - virtual ~LTI_StokesEinstein( ) { } - - 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them - */ - void getMatrixTransProp( DenseMatrix &mat, doublereal* speciesValues = 0 ) ; - protected: - - std::vector m_viscosity; - std::vector m_hydroRadius; - - }; - - //! 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 ) ; - - //! Return the matrix of binary interaction parameters. - /** - * Takes the proper mixing rule for the binary interaction parameters - * and calculates them: Not Implemented for this mixing rule - */ - 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/Makefile.in b/Cantera/src/transport/Makefile.in index 6fa9b4734..e21d70d01 100644 --- a/Cantera/src/transport/Makefile.in +++ b/Cantera/src/transport/Makefile.in @@ -36,12 +36,12 @@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG) # Base Transport Object Files TRAN_OBJ = TransportFactory.o MultiTransport.o MixTransport.o MMCollisionInt.o \ SolidTransport.o DustyGasTransport.o TransportBase.o WaterTransport.o \ - SimpleTransport.o LiquidTransportData.o LiquidTransportParams.o + SimpleTransport.o LiquidTransportData.o LiquidTransportParams.o LiquidTranInteraction.o TRAN_H = TransportFactory.h MultiTransport.h MixTransport.h \ MMCollisionInt.h SolidTransport.h DustyGasTransport.h \ TransportBase.h L_matrix.h TransportParams.h WaterTransport.h \ - SimpleTransport.h + SimpleTransport.h LiquidTranInteraction.h ifeq ($(do_electro),1) do_issp = 1 diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index aed521510..3ae2482ee 100644 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -6,6 +6,10 @@ * * */ +/* + * $Revision$ + * $Date$ + */ // turn off warnings under Windows #ifdef WIN32 @@ -37,6 +41,7 @@ #include "XML_Writer.h" #include "TransportParams.h" #include "LiquidTransportParams.h" +#include "LiquidTranInteraction.h" #include "global.h" #include "IdealGasPhase.h" #include "ctml.h" @@ -342,54 +347,54 @@ namespace Cantera { transport properties are addressed by the LTPspecies returned by newLTP */ - LiquidTranInteraction* TransportFactory::newLTI( const XML_Node &trNode, - TransportPropertyList tp_ind, - LiquidTransportParams& trParam) { + LiquidTranInteraction* TransportFactory::newLTI(const XML_Node &trNode, + TransportPropertyList tp_ind, + LiquidTransportParams& trParam) { LiquidTranInteraction* lti = 0; thermo_t* thermo = trParam.thermo; std::string model = trNode["model"]; - switch ( m_LTImodelMap[model] ) { + switch (m_LTImodelMap[model] ) { case LTI_MODEL_SOLVENT: - lti = new LTI_Solvent( tp_ind ); - lti->init( trNode, thermo ); + lti = new LTI_Solvent(tp_ind); + lti->init(trNode, thermo ); break; case LTI_MODEL_MOLEFRACS: - lti = new LTI_MoleFracs( tp_ind ); - lti->init( trNode, thermo ); + lti = new LTI_MoleFracs(tp_ind ); + lti->init(trNode, thermo ); break; case LTI_MODEL_MASSFRACS: - lti = new LTI_MassFracs( tp_ind ); - lti->init( trNode, thermo ); + lti = new LTI_MassFracs(tp_ind ); + lti->init(trNode, thermo ); break; case LTI_MODEL_LOG_MOLEFRACS: - lti = new LTI_Log_MoleFracs( tp_ind ); - lti->init( trNode, thermo ); + lti = new LTI_Log_MoleFracs(tp_ind ); + lti->init(trNode, thermo ); break; case LTI_MODEL_PAIRWISE_INTERACTION: - lti = new LTI_Pairwise_Interaction( tp_ind ); - lti->init( trNode, thermo ); - lti->setParameters( trParam ); + lti = new LTI_Pairwise_Interaction(tp_ind ); + 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 ); + 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 ); + 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 ); + 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 ); - lti->init( trNode, thermo ); + lti = new LiquidTranInteraction(tp_ind ); + lti->init(trNode, thermo ); } return lti; } @@ -1002,7 +1007,7 @@ namespace Cantera { switch (m_tranPropMap[nodeName]) { case TP_VISCOSITY: - data.viscosity = newLTP(xmlChild, name, m_tranPropMap[nodeName], temp_thermo ); + data.viscosity = newLTP(xmlChild, name, m_tranPropMap[nodeName], temp_thermo); break; case TP_IONCONDUCTIVITY: data.ionConductivity = newLTP(xmlChild, name, m_tranPropMap[nodeName], temp_thermo); @@ -1034,25 +1039,25 @@ namespace Cantera { data.thermalCond = newLTP(xmlChild, name, m_tranPropMap[nodeName], - temp_thermo ); + temp_thermo); break; case TP_DIFFUSIVITY: data.speciesDiffusivity = newLTP(xmlChild, name, m_tranPropMap[nodeName], - temp_thermo ); + temp_thermo); break; case TP_HYDRORADIUS: data.hydroRadius = newLTP(xmlChild, name, m_tranPropMap[nodeName], - temp_thermo ); + temp_thermo); break; case TP_ELECTCOND: data.electCond = newLTP(xmlChild, name, m_tranPropMap[nodeName], - temp_thermo ); + temp_thermo); break; default: @@ -1126,9 +1131,9 @@ namespace Cantera { trParam.viscosity = newLTI(compDepNode, m_tranPropMap[nodeName], trParam); break; case TP_IONCONDUCTIVITY: - trParam.ionConductivity = newLTI( compDepNode, + trParam.ionConductivity = newLTI(compDepNode, m_tranPropMap[nodeName], - trParam ); + trParam); break; case TP_MOBILITYRATIO: { @@ -1139,9 +1144,9 @@ namespace Cantera { string firstSpec = specName.substr(0,loc); string secondSpec = specName.substr(loc+1); int index = temp_thermo->speciesIndex(firstSpec.c_str())+nsp*temp_thermo->speciesIndex(secondSpec.c_str()); - trParam.mobilityRatio[index] = newLTI( propSpecNode, + trParam.mobilityRatio[index] = newLTI(propSpecNode, m_tranPropMap[nodeName], - trParam ); + trParam); }; }; break; @@ -1151,21 +1156,21 @@ namespace Cantera { XML_Node &propSpecNode = compDepNode.child(iSpec); string specName = propSpecNode.name(); int index = temp_thermo->speciesIndex(specName.c_str()); - trParam.selfDiffusion[index] = newLTI( propSpecNode, + trParam.selfDiffusion[index] = newLTI(propSpecNode, m_tranPropMap[nodeName], - trParam ); + trParam); }; }; break; case TP_THERMALCOND: - trParam.thermalCond = newLTI( compDepNode, + trParam.thermalCond = newLTI(compDepNode, m_tranPropMap[nodeName], - trParam ); + trParam); break; case TP_DIFFUSIVITY: - trParam.speciesDiffusivity = newLTI( compDepNode, + trParam.speciesDiffusivity = newLTI(compDepNode, m_tranPropMap[nodeName], - trParam ); + trParam); break; case TP_HYDRORADIUS: trParam.hydroRadius = newLTI(compDepNode, @@ -1173,12 +1178,12 @@ namespace Cantera { trParam); break; case TP_ELECTCOND: - trParam.electCond = newLTI( compDepNode, + trParam.electCond = newLTI(compDepNode, m_tranPropMap[nodeName], - trParam ); + trParam); break; default: - throw CanteraError("getLiquidInteractionsTransportData","unknown transport property: " + nodeName ); + throw CanteraError("getLiquidInteractionsTransportData","unknown transport property: " + nodeName); } } @@ -1201,7 +1206,7 @@ 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. "); } } } @@ -1320,7 +1325,7 @@ namespace Cantera { // self-diffusion coefficient, without polar // corrections diffcoeff = ThreeSixteenths * - sqrt( 2.0 * Pi/tr.reducedMass(k,k) ) * + sqrt(2.0 * Pi/tr.reducedMass(k,k)) * pow((Boltzmann * t), 1.5)/ (Pi * tr.sigma[k] * tr.sigma[k] * om11); @@ -1483,7 +1488,7 @@ namespace Cantera { om11 = m_integrals->omega11(tstar, tr.delta(j,k)); diffcoeff = ThreeSixteenths * - sqrt( 2.0 * Pi/tr.reducedMass(k,j) ) * + sqrt(2.0 * Pi/tr.reducedMass(k,j)) * pow((Boltzmann * t), 1.5)/ (Pi * sigma * sigma * om11); diff --git a/Cantera/src/transport/TransportFactory.h b/Cantera/src/transport/TransportFactory.h index fdb55db59..504c727c9 100644 --- a/Cantera/src/transport/TransportFactory.h +++ b/Cantera/src/transport/TransportFactory.h @@ -12,8 +12,8 @@ * */ -#ifndef CT_TRANFACTORY_H -#define CT_TRANFACTORY_H +#ifndef CT_TRANSPORTFACTORY_H +#define CT_TRANSPORTFACTORY_H // turn off warnings under Windows @@ -137,17 +137,20 @@ namespace Cantera { newLTP( const XML_Node &trNode, std::string &name, TransportPropertyList tp_ind, thermo_t* thermo) ; - - /** - * make one of several transport models, and return a base class - * pointer to it. This method operates at the level of a - * single mixture transport property. Individual species - * transport properties are addressed by the LTPspecies - * returned by newLTP + + //! Factory function for the construction of new LiquidTranInteraction + //! objects, which are transport models. + /*! + * This method operates at the level of a single mixture transport property. Individual species + * transport properties are addressed by the LTPspecies returned by newLTP. + * + * @param trNode XML_Node containing the information for the interaction + * @param tp_ind TransportPropertylist object + * @param trParam reference to the LiquidTransportParams object */ - virtual LiquidTranInteraction* newLTI( const XML_Node &trNode, - TransportPropertyList tp_ind, - LiquidTransportParams& trParam) ; + virtual LiquidTranInteraction* newLTI(const XML_Node &trNode, + TransportPropertyList tp_ind, + LiquidTransportParams& trParam); //! Build a new transport manager using a transport manager diff --git a/Cantera/src/transport/TransportParams.h b/Cantera/src/transport/TransportParams.h index de03f32bd..9eb85bd9c 100644 --- a/Cantera/src/transport/TransportParams.h +++ b/Cantera/src/transport/TransportParams.h @@ -11,6 +11,13 @@ namespace Cantera { + class NotImplemented : public CanteraError { + public: + NotImplemented(std::string method) : CanteraError("Transport", + "\n\n**** Method "+method+" not implemented. ****\n" + "(Did you forget to specify a transport model?)\n\n") {} + }; + /** * Base class to hold transport model parameters. * Used by TransportFactory.