diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index 8175c32eb..08b566dea 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -443,35 +443,6 @@ namespace Cantera { } } - void LiquidTransport::getSpeciesDiffusiveMassFluxes(doublereal* const fluxes) { - int n, k; - - update_temp(); - update_conc(); - - - getMixDiffCoeffs(DATA_PTR(m_spwork)); - - - const array_fp& mw = m_thermo->molecularWeights(); - const doublereal* const y = m_thermo->massFractions(); - const doublereal rhon = m_thermo->molarDensity(); - // Unroll wrt ndim - vector_fp sum(m_nDim,0.0); - for (n = 0; n < m_nDim; n++) { - for (k = 0; k < m_nsp; k++) { - fluxes[n*m_nsp + k] = -rhon * mw[k] * m_spwork[k] * m_Grad_X[n*m_nsp + k]; - sum[n] += fluxes[n*m_nsp + k]; - } - } - // add correction flux to enforce sum to zero - for (n = 0; n < m_nDim; n++) { - for (k = 0; k < m_nsp; k++) { - fluxes[n*m_nsp + k] -= y[k]*sum[n]; - } - } - } - /** * Mixture-averaged diffusion coefficients [m^2/s]. * diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h index 852641192..da8a1751b 100644 --- a/Cantera/src/transport/LiquidTransport.h +++ b/Cantera/src/transport/LiquidTransport.h @@ -326,19 +326,6 @@ namespace Cantera { int ldx, const doublereal* grad_X, int ldf, doublereal* fluxes); - //! Return the species diffusive mass fluxes - /*! - * - * - * - * @param ndim The number of spatial dimensions (1, 2, or 3). - * @param grad_T The temperature gradient (ignored in this model). - * @param ldx Leading dimension of the grad_X array. - * The diffusive mass flux of species \e k is computed from - * - * - */ - virtual void getSpeciesDiffusiveMassFluxes(doublereal* const fluxes); /** * @param ndim The number of spatial dimensions (1, 2, or 3). @@ -719,14 +706,15 @@ namespace Cantera { int m_nDim; private: - - /** - * Throw an exception if this method is invoked. + + //! Throw an exception if this method is invoked. + /*! * This probably indicates something is not yet implemented. + * + * @pram msg Indicates the member function which is not implemented */ doublereal err(std::string msg) const; - }; } #endif diff --git a/Cantera/src/transport/LiquidTransportData.h b/Cantera/src/transport/LiquidTransportData.h new file mode 100644 index 000000000..f5b1bd248 --- /dev/null +++ b/Cantera/src/transport/LiquidTransportData.h @@ -0,0 +1,80 @@ +/** + * @file TransportFactory.h + * Header file defining class TransportFactory + * (see \link Cantera::TransportFactory TransportFactory\endlink) + */ +/* + * $Author: hkmoffa $ + * $Date: 2008/12/24 18:19:01 $ + * $Revision: 1.14 $ + * + * Copyright 2001 California Institute of Technology + * + */ + +#ifndef CT_LIQUIDTRANSPORTDATA_H +#define CT_LIQUIDTRANSPORTDATA_H + + +// STL includes +#include +#include +#include +#include + + + +// Cantera includes +#include "ct_defs.h" +#include "TransportBase.h" +#include "FactoryBase.h" + + +namespace Cantera { + + enum LiquidTR_Model { + LTR_MODEL_NOTSET=-1, + LTR_MODEL_CONSTANT, + LTR_MODEL_ARRHENIUS, + LTR_MODEL_COEFF + }; + + class LiquidTransportData { + + public: + + LiquidTransportData() : + speciesName("-"), + model_hydroradius(LTR_MODEL_NOTSET), + hydroradius(-1.0), + model_viscosity(LTR_MODEL_NOTSET), + model_thermalCond(LTR_MODEL_NOTSET), + model_speciesDiffusivity(LTR_MODEL_NOTSET) + { + } + + std::string speciesName; + + //! Model type for the hydroradius + LiquidTR_Model model_hydroradius; + + //! Actual value of the hydroradius + doublereal hydroradius; + + //! Model type for the hydroradius + LiquidTR_Model model_viscosity; + vector_fp viscCoeffs; + + //! Model type for the hydroradius + LiquidTR_Model model_thermalCond; + + vector_fp thermalCondCoeffs; + + //! Model type for the hydroradius + LiquidTR_Model model_speciesDiffusivity; + + vector_fp speciesDiffusivityCoeffs; + }; + +} +#endif diff --git a/Cantera/src/transport/LiquidTransportParams.h b/Cantera/src/transport/LiquidTransportParams.h index 08142b745..9b0d98c94 100644 --- a/Cantera/src/transport/LiquidTransportParams.h +++ b/Cantera/src/transport/LiquidTransportParams.h @@ -6,6 +6,7 @@ #include "ct_defs.h" #include "TransportBase.h" #include "TransportParams.h" +#include "LiquidTransportData.h" #include "xml.h" #include "XML_Writer.h" @@ -81,6 +82,8 @@ namespace Cantera { vector_fp B_k_cond; + std::vector LTData; + }; } diff --git a/Cantera/src/transport/Makefile.in b/Cantera/src/transport/Makefile.in index ad6d76b15..9f52f08e1 100644 --- a/Cantera/src/transport/Makefile.in +++ b/Cantera/src/transport/Makefile.in @@ -35,11 +35,13 @@ 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 + SolidTransport.o DustyGasTransport.o TransportBase.o WaterTransport.o \ + SimpleTransport.o TRAN_H = TransportFactory.h MultiTransport.h MixTransport.h \ MMCollisionInt.h SolidTransport.h DustyGasTransport.h \ - TransportBase.h L_matrix.h TransportParams.h WaterTransport.h + TransportBase.h L_matrix.h TransportParams.h WaterTransport.h \ + SimpleTransport.h LiquidTransportData.h ifeq ($(do_electro),1) do_issp = 1 diff --git a/Cantera/src/transport/MixTransport.h b/Cantera/src/transport/MixTransport.h old mode 100755 new mode 100644 diff --git a/Cantera/src/transport/MultiTransport.h b/Cantera/src/transport/MultiTransport.h old mode 100755 new mode 100644 diff --git a/Cantera/src/transport/SimpleTransport.cpp b/Cantera/src/transport/SimpleTransport.cpp new file mode 100644 index 000000000..ed119bfe5 --- /dev/null +++ b/Cantera/src/transport/SimpleTransport.cpp @@ -0,0 +1,569 @@ +/** + * @file SimpleTransport.cpp + * Simple mostly constant transport properties + */ +/* + * $Revision: 1.10 $ + * $Date: 2009/03/24 20:44:30 $ + */ + +#include "ThermoPhase.h" +#include "SimpleTransport.h" + +#include "utilities.h" +#include "LiquidTransportParams.h" +#include "TransportFactory.h" + +#include "ctlapack.h" + +#include +using namespace std; + +/** + * Mole fractions below MIN_X will be set to MIN_X when computing + * transport properties. + */ +#define MIN_X 1.e-14 + + +namespace Cantera { + //================================================================================================ + SimpleTransport::SimpleTransport(thermo_t* thermo, int ndim) : + Transport(thermo, ndim), + m_nsp(0), + m_tmin(-1.0), + m_tmax(100000.), + m_iStateMF(-1), + m_temp(-1.0), + m_press(-1.0), + m_lambda(-1.0), + m_viscmix(-1.0), + m_visc_mix_ok(false), + m_visc_temp_ok(false), + m_diff_mix_ok(false), + m_diff_temp_ok(false), + m_cond_temp_ok(false), + m_cond_mix_ok(false) + { + } + //================================================================================================ + SimpleTransport::SimpleTransport(const SimpleTransport &right) : + Transport(), + m_nsp(0), + m_tmin(-1.0), + m_tmax(100000.), + m_iStateMF(-1), + m_temp(-1.0), + m_press(-1.0), + m_lambda(-1.0), + m_viscmix(-1.0), + m_visc_mix_ok(false), + m_visc_temp_ok(false), + m_diff_mix_ok(false), + m_diff_temp_ok(false), + m_cond_temp_ok(false), + m_cond_mix_ok(false) + { + /* + * Use the assignment operator to do the brunt + * of the work for the copy construtor. + */ + *this = right; + } + //================================================================================================ + SimpleTransport& SimpleTransport::operator=(const SimpleTransport& right) { + if (&right != this) { + return *this; + } + Transport::operator=(right); + m_nsp = right.m_nsp; + m_tmin = right.m_tmin; + m_tmax = right.m_tmax; + m_mw = right.m_mw; + + m_Grad_X = right.m_Grad_X; + m_Grad_T = right.m_Grad_T; + m_Grad_V = right.m_Grad_V; + + m_viscSpecies = right.m_viscSpecies; + m_condSpecies = right.m_condSpecies; + m_iStateMF = -1; + m_molefracs = right.m_molefracs; + m_concentrations = right.m_concentrations; + m_chargeSpecies = right.m_chargeSpecies; + + + m_temp = right.m_temp; + m_press = right.m_press; + m_lambda = right.m_lambda; + m_viscmix = right.m_viscmix; + m_spwork = right.m_spwork; + m_visc_mix_ok = false; + m_visc_temp_ok = false; + m_diff_mix_ok = false; + m_diff_temp_ok = false; + m_cond_temp_ok = false; + m_cond_mix_ok = false; + m_nDim = right.m_nDim; + + return *this; + } + + //================================================================================================ + Transport *SimpleTransport::duplMyselfAsTransport() const { + SimpleTransport* tr = new SimpleTransport(*this); + return (dynamic_cast(tr)); + } + //================================================================================================ + // Initialize the object + /* + * This is where we dimension everything. + */ + bool SimpleTransport::initLiquid(LiquidTransportParams& tr) { + + // constant substance attributes + m_thermo = tr.thermo; + m_nsp = m_thermo->nSpecies(); + m_tmin = m_thermo->minTemp(); + m_tmax = m_thermo->maxTemp(); + + // make a local copy of the molecular weights + m_mw.resize(m_nsp); + copy(m_thermo->molecularWeights().begin(), + m_thermo->molecularWeights().end(), m_mw.begin()); + + //save logarithm of pre-exponential for easier computation + + //m_diffcoeffs = tr.diffcoeffs; + + + m_viscSpecies.resize(m_nsp); + m_condSpecies.resize(m_nsp); + + + m_molefracs.resize(m_nsp); + m_spwork.resize(m_nsp); + + // resize the internal gradient variables + m_Grad_X.resize(m_nDim * m_nsp, 0.0); + m_Grad_T.resize(m_nDim, 0.0); + m_Grad_V.resize(m_nDim, 0.0); + + + + // set all flags to false + m_visc_mix_ok = false; + m_visc_temp_ok = false; + + m_cond_temp_ok = false; + m_cond_mix_ok = false; + + m_diff_temp_ok = false; + m_diff_mix_ok = false; + + return true; + } + + //================================================================================================ + // Returns the mixture viscosity of the solution + /* + * The viscosity is computed using the general mixture rules + * specified in the variable compositionDepType_. + * + * Solvent-only: + * \f[ + * \mu = \mu_0 + * \f] + * Mixture-average: + * \f[ + * \mu = \sum_k {\mu_k X_k} + * \f] + * + * Here \f$ \mu_k \f$ is the viscosity of pure species \e k. + * + * @see updateViscosity_T(); + */ + doublereal SimpleTransport::viscosity() { + + update_T(); + update_C(); + + if (m_visc_mix_ok) return m_viscmix; + + // update m_viscSpecies[] if necessary + if (!m_visc_temp_ok) { + updateViscosity_T(); + } + + if (compositionDepType_ == 0) { + m_viscmix = m_viscSpecies[0]; + } else if (compositionDepType_ == 1) { + m_viscmix = 0.0; + for (int k = 0; k < m_nsp; k++) { + m_viscmix += m_viscSpecies[k] * m_molefracs[k]; + } + } + m_visc_mix_ok = true; + return m_viscmix; + } + //================================================================================================ + void SimpleTransport::getSpeciesViscosities(doublereal* visc) { + update_T(); + if (!m_visc_temp_ok) { + updateViscosity_T(); + } + copy(m_viscSpecies.begin(), m_viscSpecies.end(), visc); + } + //================================================================================================ + void SimpleTransport::getBinaryDiffCoeffs(int ld, doublereal* d) { + int i, j; + double bdiff; + update_T(); + + // if necessary, evaluate the species diffusion coefficents + // from the polynomial fits + if (!m_diff_temp_ok) updateDiff_T(); + + for (i = 0; i < m_nsp; i++) { + for (j = 0; j < m_nsp; j++) { + bdiff = 0.5 * (m_diffSpecies[i] + m_diffSpecies[j]); + d[i*m_nsp+j] = bdiff; + } + } + } + //================================================================================================ + void SimpleTransport::getMobilities(doublereal* const mobil) { + // this needs to be checked out. + int k; + getMixDiffCoeffs(DATA_PTR(m_spwork)); + doublereal c1 = ElectronCharge / (Boltzmann * m_temp); + for (k = 0; k < m_nsp; k++) { + mobil[k] = c1 * m_spwork[k] * m_thermo->charge(k); + } + } + //================================================================================================ + void SimpleTransport::set_Grad_V(const doublereal* const grad_V) { + for (int a = 0; a < m_nDim; a++) { + m_Grad_V[a] = grad_V[a]; + } + } + //================================================================================================ + void SimpleTransport::set_Grad_T(const doublereal* const grad_T) { + for (int a = 0; a < m_nDim; a++) { + m_Grad_T[a] = grad_T[a]; + } + } + //================================================================================================ + void SimpleTransport::set_Grad_X(const doublereal* const grad_X) { + int itop = m_nDim * m_nsp; + for (int i = 0; i < itop; i++) { + m_Grad_X[i] = grad_X[i]; + } + } + + //================================================================================================ + // Returns the mixture thermal conductivity of the solution + /* + * The thermal is computed using the general mixture rules + * specified in the variable compositionDepType_. + * + * Solvent-only: + * \f[ + * \lambda = \lambda_0 + * \f] + * Mixture-average: + * \f[ + * \lambda = \sum_k {\lambda_k X_k} + * \f] + * + * Here \f$ \lambda_k \f$ is the thermal conductivity of pure species \e k. + * + * @see updateCond_T(); + */ + doublereal SimpleTransport::thermalConductivity() { + update_T(); + update_C(); + if (!m_cond_temp_ok) { + updateCond_T(); + } + if (!m_cond_mix_ok) { + if (compositionDepType_ == 0) { + m_lambda = m_condSpecies[0]; + } else if (compositionDepType_ == 1) { + m_lambda = 0.0; + for (int k = 0; k < m_nsp; k++) { + m_lambda += m_condSpecies[k] * m_molefracs[k]; + } + } + m_cond_mix_ok = true; + } + return m_lambda; + } + //================================================================================================ + + /* + * Thermal diffusion is not considered in this mixture-averaged + * model. To include thermal diffusion, use transport manager + * MultiTransport instead. This methods fills out array dt with + * zeros. + */ + void SimpleTransport::getThermalDiffCoeffs(doublereal* const dt) { + for (int k = 0; k < m_nsp; k++) { + dt[k] = 0.0; + } + } +//================================================================================================ + /** + * @param ndim The number of spatial dimensions (1, 2, or 3). + * @param grad_T The temperature gradient (ignored in this model). + * @param ldx Leading dimension of the grad_X array. + * The diffusive mass flux of species \e k is computed from + * + * \f[ + * \vec{j}_k = -n M_k D_k \nabla X_k. + * \f] + */ + void SimpleTransport::getSpeciesFluxes(int ndim, + const doublereal* grad_T, + int ldx, const doublereal* grad_X, + int ldf, doublereal* fluxes) { + set_Grad_T(grad_T); + set_Grad_X(grad_X); + getSpeciesFluxesExt(ldf, fluxes); + } + //================================================================================================ + // Return the species diffusive mass fluxes wrt to + // the mass averaged velocity, + /* + * + * units = kg/m2/s + * + * Internally, gradients in the in mole fraction, temperature + * and electrostatic potential contribute to the diffusive flux + * + * + * The diffusive mass flux of species \e k is computed from the following + * formula + * + * \f[ + * j_k = - \rho M_k D_k \nabla X_k - Y_k V_c + * \f] + * + * where V_c is the correction velocity + * + * \f[ + * V_c = - \sum_j {\rho M_j D_j \nabla X_j} + * \f] + * + * @param ldf stride of the fluxes array. Must be equal to + * or greater than the number of species. + * @param fluxes Vector of calculated fluxes + */ + void SimpleTransport::getSpeciesFluxesExt(int ldf, doublereal* fluxes) { + int n, k; + AssertThrow(ldf >= m_nsp ,"SimpleTransport::getSpeciesFluxesExt: Stride must be greater than m_nsp"); + update_T(); + update_C(); + + getMixDiffCoeffs(DATA_PTR(m_spwork)); + + + const array_fp& mw = m_thermo->molecularWeights(); + const doublereal* y = m_thermo->massFractions(); + doublereal rhon = m_thermo->molarDensity(); + // Unroll wrt ndim + vector_fp sum(m_nDim,0.0); + for (n = 0; n < m_nDim; n++) { + for (k = 0; k < m_nsp; k++) { + fluxes[n*ldf + k] = -rhon * mw[k] * m_spwork[k] * m_Grad_X[n*m_nsp + k]; + sum[n] += fluxes[n*ldf + k]; + } + } + // add correction flux to enforce sum to zero + for (n = 0; n < m_nDim; n++) { + for (k = 0; k < m_nsp; k++) { + fluxes[n*ldf + k] -= y[k]*sum[n]; + } + } + } + //================================================================================================ + // Mixture-averaged diffusion coefficients [m^2/s]. + /* + * Returns the simple diffusion coefficients input into the model. Nothing fancy here. + */ + void SimpleTransport::getMixDiffCoeffs(doublereal* const d) { + update_T(); + update_C(); + // update the binary diffusion coefficients if necessary + if (!m_diff_temp_ok) { + updateDiff_T(); + } + for (int k = 0; k < m_nsp; k++) { + d[k] = m_diffSpecies[k]; + } + } + //================================================================================================ + + // Handles the effects of changes in the mixture concentration + /* + * This is called for every interface call to check whether + * the concentrations have changed. Concentrations change + * whenever the pressure or the mole fraction has changed. + * If it has changed, the recalculations should be done. + * + * Note this should be a lightweight function since it's + * part of all of the interfaces. + * + * @internal + */ + bool SimpleTransport::update_C() { + // If the pressure has changed then the concentrations + // have changed. + doublereal pres = m_thermo->pressure(); + bool qReturn = true; + if (pres != m_press) { + qReturn = false; + m_press = pres; + } + int iStateNew = m_thermo->stateMFNumber(); + if (iStateNew != m_iStateMF) { + qReturn = false; + m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); + m_thermo->getConcentrations(DATA_PTR(m_concentrations)); + concTot_ = 0.0; + for (int k = 0; k < m_nsp; k++) { + m_molefracs[k] = fmaxx(0.0, m_molefracs[k]); + concTot_ += m_concentrations[k]; + } + dens_ = m_thermo->density(); + meanMolecularWeight_ = m_thermo->meanMolecularWeight(); + } + if (qReturn) { + return false; + } + + + // Mixture stuff needs to be evaluated + m_visc_mix_ok = false; + m_diff_mix_ok = false; + m_cond_mix_ok = false; + + return true; + } + + //================================================================================================ + /** + * Update the temperature-dependent parts of the mixture-averaged + * thermal conductivity. + */ + void SimpleTransport::updateCond_T() { + int k; + if (tempDepType_ == 0) { + for (k = 0; k < m_nsp; k++) { + Coeff_T_ &coeff = m_coeffLambda_Ns[k]; + m_condSpecies[k] = coeff[0]; + } + } else if (tempDepType_ == 1) { + for (k = 0; k < m_nsp; k++) { + Coeff_T_ &coeff = m_coeffLambda_Ns[k]; + m_condSpecies[k] = coeff[0] * pow(m_temp,coeff[1]) * exp(-coeff[2]/m_temp); + } + } + m_cond_temp_ok = true; + m_cond_mix_ok = false; + } + //================================================================================================ + /** + * Update the species diffusion coefficients. + */ + void SimpleTransport::updateDiff_T() { + int k; + if (tempDepType_ == 0) { + for (k = 0; k < m_nsp; k++) { + Coeff_T_ &coeff = m_coeffDiff_Ns[k]; + m_diffSpecies[k] = coeff[0]; + } + } else if (tempDepType_ == 1) { + for (k = 0; k < m_nsp; k++) { + Coeff_T_ &coeff = m_coeffDiff_Ns[k]; + m_viscSpecies[k] = coeff[0] * pow(m_temp,coeff[1]) * exp(-coeff[2]/m_temp); + } + } + m_diff_temp_ok = true; + m_diff_mix_ok = false; + } + //================================================================================================ + + /** + * Update the pure-species viscosities. + */ + void SimpleTransport::updateViscosities_C() { + + } + //================================================================================================ + /** + * Update the temperature-dependent viscosity terms. + * Updates the array of pure species viscosities, and the + * weighting functions in the viscosity mixture rule. + * The flag m_visc_ok is set to true. + */ + void SimpleTransport::updateViscosity_T() { + int k; + if (tempDepType_ == 0) { + for (k = 0; k < m_nsp; k++) { + Coeff_T_ &coeff = m_coeffVisc_Ns[k]; + m_viscSpecies[k] = coeff[0]; + } + } else if (tempDepType_ == 1) { + for (k = 0; k < m_nsp; k++) { + Coeff_T_ &coeff = m_coeffVisc_Ns[k]; + m_viscSpecies[k] = coeff[0] * pow(m_temp,coeff[1]) * exp(-coeff[2]/m_temp); + } + } + m_visc_temp_ok = true; + m_visc_mix_ok = false; + } + //================================================================================================= + bool SimpleTransport::update_T() + { + doublereal t = m_thermo->temperature(); + if (t == m_temp) return false; + if (t < 0.0) { + throw CanteraError("SimpleTransport::update_T", + "negative temperature "+fp2str(t)); + } + + // Compute various functions of temperature + m_temp = t; + + // temperature has changed, so polynomial temperature + // interpolations will need to be reevaluated. + // Set all of these flags to false + m_visc_mix_ok = false; + m_visc_temp_ok = false; + + m_cond_temp_ok = true; + m_cond_mix_ok = false; + + m_diff_mix_ok = false; + m_diff_temp_ok = false; + + return true; + } + + /** + * Throw an exception if this method is invoked. + * This probably indicates something is not yet implemented. + */ + doublereal SimpleTransport::err(std::string msg) const { + throw CanteraError("SimpleTransport Class", + "\n\n\n**** Method "+ msg +" not implemented in model " + + int2str(model()) + " ****\n" + "(Did you forget to specify a transport model?)\n\n\n"); + + return 0.0; + } + //================================================================================================ + +} +//================================================================================================ diff --git a/Cantera/src/transport/SimpleTransport.h b/Cantera/src/transport/SimpleTransport.h new file mode 100644 index 000000000..45f828564 --- /dev/null +++ b/Cantera/src/transport/SimpleTransport.h @@ -0,0 +1,666 @@ +/** + * + * @file SimpleTransport.h + * Header file defining class SimpleTransport + */ +/* + * $Revision: 1.9 $ + * $Date: 2009/03/27 18:24:39 $ + */ + +#ifndef CT_SIMPLETRAN_H +#define CT_SIMPLETRAN_H + + + +// STL includes +#include +#include +#include +#include +#include + +using namespace std; + +// Cantera includes +#include "TransportBase.h" +#include "DenseMatrix.h" +#include "TransportParams.h" +#include "LiquidTransportParams.h" + +namespace Cantera { + + + + + class LiquidTransportParams; + + + //! Class LiquidTransport implements mixture-averaged transport + //! properties for liquid phases. + /*! + * The model is based on that + * described by Newman, Electrochemical Systems + * + * The velocity of species i may be described by the + * following equation p. 297 (12.1) + * + * \f[ + * c_i \nabla \mu_i = R T \sum_j \frac{c_i c_j}{c_T D_{ij}} + * (\mathbf{v}_j - \mathbf{v}_i) + * \f] + * + * This as written is degenerate by 1 dof. + * + * To fix this we must add in the definition of the mass averaged + * velocity of the solution. We will call the simple bold-faced + * \f$\mathbf{v} \f$ + * symbol the mass-averaged velocity. Then, the relation + * between \f$\mathbf{v}\f$ and the individual species velocities is + * \f$\mathbf{v}_i\f$ + * + * \f[ + * \rho_i \mathbf{v}_i = \rho_i \mathbf{v} + \mathbf{j}_i + * \f] + * where \f$\mathbf{j}_i\f$ are the diffusional fluxes of species i + * with respect to the mass averaged velocity and + * + * \f[ + * \sum_i \mathbf{j}_i = 0 + * \f] + * + * and + * + * \f[ + * \sum_i \rho_i \mathbf{v}_i = \rho \mathbf{v} + * \f] + * + * Using these definitions, we can write + * + * \f[ + * \mathbf{v}_i = \mathbf{v} + \frac{\mathbf{j}_i}{\rho_i} + * \f] + * + * + * \f[ + * c_i \nabla \mu_i = R T \sum_j \frac{c_i c_j}{c_T D_{ij}} + * (\frac{\mathbf{j}_j}{\rho_j} - \frac{\mathbf{j}_i}{\rho_i}) + * = R T \sum_j \frac{1}{D_{ij}} + * (\frac{x_i \mathbf{j}_j}{M_j} - \frac{x_j \mathbf{j}_i}{M_i}) + * \f] + * + * The equations that we actually solve are + * + * \f[ + * c_i \nabla \mu_i = + * = R T \sum_j \frac{1}{D_{ij}} + * (\frac{x_i \mathbf{j}_j}{M_j} - \frac{x_j \mathbf{j}_i}{M_i}) + * \f] + * and we replace the 0th equation with the following: + * + * \f[ + * \sum_i \mathbf{j}_i = 0 + * \f] + * + * When there are charged species, we replace the rhs with the + * gradient of the electrochemical potential to obtain the + * modified equation + * + * \f[ + * c_i \nabla \mu_i + c_i F z_i \nabla \Phi + * = R T \sum_j \frac{1}{D_{ij}} + * (\frac{x_i \mathbf{j}_j}{M_j} - \frac{x_j \mathbf{j}_i}{M_i}) + * \f] + * + * With this formulation we may solve for the diffusion velocities, + * without having to worry about what the mass averaged velocity + * is. + * + *

Viscosity Calculation

+ * + * The viscosity calculation may be broken down into two parts. + * In the first part, the viscosity of the pure species are calculated + * In the second part, a mixing rule is applied, based on the + * Wilkes correlation, to yield the mixture viscosity. + * + * + * + */ + class SimpleTransport : public Transport { + public: + + typedef double Coeff_T_ [4]; + + + //! Default constructor. + /*! + * This requires call to initLiquid(LiquidTransportParams& tr) + * after filling LiquidTransportParams to complete instantiation. + * The filling of LiquidTransportParams is currently carried out + * in the TransportFactory class, but might be moved at some point. + * + * @param thermo ThermoPhase object holding species information. + * @param ndim Number of spatial dimensions. + */ + SimpleTransport(thermo_t* thermo = 0, int ndim = 1); + + //!Copy Constructor for the %LiquidThermo object. + /*! + * @param right %LiquidTransport to be copied + */ + SimpleTransport(const SimpleTransport &right); + + //! Assignment operator + /*! + * This is NOT a virtual function. + * + * @param right Reference to %LiquidTransport object to be copied + * into the current one. + */ + SimpleTransport& operator=(const SimpleTransport& right); + + //! Duplication routine for objects which inherit from + //! %Transport + /*! + * This virtual routine can be used to duplicate %Transport objects + * inherited from %Transport even if the application only has + * a pointer to %Transport to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + */ + virtual Transport *duplMyselfAsTransport() const; + + + //! virtual destructor + virtual ~SimpleTransport() {} + + //! Initialize the transport object + /*! + * Here we change all of the internal dimensions to be sufficient. + * We get the object ready to do property evaluations. + * + * @param tr Transport parameters for all of the species + * in the phase. + */ + virtual bool initLiquid(LiquidTransportParams& tr); + + friend class TransportFactory; + + + //! Return the model id for this transport parameterization + virtual int model() const { + return cSimpleTransport; + } + + //! overloaded base class methods + + //! Returns the mixture viscosity of the solution + /*! + * The viscosity is computed using the general mixture rules + * specified in the variable compositionDepType_. + * + * Solvent-only: + * \f[ + * \mu = \mu_0 + * \f] + * Mixture-average: + * \f[ + * \mu = \sum_k {\mu_k X_k} + * \f] + * + * Here \f$ \mu_k \f$ is the viscosity of pure species \e k. + * + * @see updateViscosity_T(); + */ + virtual doublereal viscosity(); + + //! Returns the pure species viscosities + /*! + * The pure species viscosities are to be given in an Arrhenius + * form in accordance with activated-jump-process dominated transport. + */ + virtual void getSpeciesViscosities(doublereal* const visc); + + //! Returns the binary diffusion coefficients + /*! + * @param ld + * @param d + */ + virtual void getBinaryDiffCoeffs(const int ld, doublereal* const d); + + //! Get the Mixture diffusion coefficients + /*! + * @param d vector of mixture diffusion coefficients + * units = m2 s-1. length = number of species + */ + virtual void getMixDiffCoeffs(doublereal* const d); + + + virtual void getThermalDiffCoeffs(doublereal* const dt); + + + //! Returns the mixture thermal conductivity of the solution + /*! + * The thermal is computed using the general mixture rules + * specified in the variable compositionDepType_. + * + * Controlling update boolean = m_condmix_ok + * + * Units are in W/m/K or equivalently kg m / s3 / K + * + * Solvent-only: + * \f[ + * \lambda = \lambda_0 + * \f] + * Mixture-average: + * \f[ + * \lambda = \sum_k {\lambda_k X_k} + * \f] + * + * Here \f$ \lambda_k \f$ is the thermal conductivity of pure species \e k. + * + * @see updateCond_T(); + */ + virtual doublereal thermalConductivity(); + + //! Get the Mobilities + /*! + * @param mobil + */ + virtual void getMobilities(doublereal* const mobil); + + //! Specify the valpdaue of the gradient of the voltage + /*! + * + * @param grad_V Gradient of the voltage (length num dimensions); + */ + virtual void set_Grad_V(const doublereal* const grad_V); + + //! Specify the value of the gradient of the temperature + /*! + * + * @param grad_V Gradient of the temperature (length num dimensions); + */ + virtual void set_Grad_T(const doublereal* const grad_T); + + //! Specify the value of the gradient of the MoleFractions + /*! + * + * @param grad_X Gradient of the mole fractions(length nsp * num dimensions); + */ + virtual void set_Grad_X(const doublereal* const grad_X); + + + /** + * @param ndim The number of spatial dimensions (1, 2, or 3). + * @param grad_T The temperature gradient (ignored in this model). + * @param ldx Leading dimension of the grad_X array. + * The diffusive mass flux of species \e k is computed from + * + * + */ + virtual void getSpeciesFluxes(int ndim, + const doublereal* grad_T, + int ldx, const doublereal* grad_X, + int ldf, doublereal* fluxes); + + //! Return the species diffusive mass fluxes wrt to + //! the mass averaged velocity, + /*! + * + * units = kg/m2/s + * + * Internally, gradients in the in mole fraction, temperature + * and electrostatic potential contribute to the diffusive flux + * + * + * The diffusive mass flux of species \e k is computed from the following + * formula + * + * \f[ + * j_k = - \rho M_k D_k \nabla X_k - Y_k V_c + * \f] + * + * where V_c is the correction velocity + * + * \f[ + * V_c = - \sum_j {\rho M_j D_j \nabla X_j} + * \f] + * + * @param ldf stride of the fluxes array. Must be equal to + * or greater than the number of species. + * @param fluxes Vector of calculated fluxes + */ + virtual void getSpeciesFluxesExt(int ldf, doublereal* fluxes); + + protected: + + //! Handles the effects of changes in the Temperature, internally + //! within the object. + /*! + * This is called whenever a transport property is requested. + * The first task is to check whether the temperature has changed + * since the last call to update_T(). + * If it hasn't then an immediate return is carried out. + * + * @internal + * + * @return Returns true if the temperature has changed, and false otherwise + */ + virtual bool update_T(); + + //! Handles the effects of changes in the mixture concentration + /*! + * This is called for every interface call to check whether + * the concentrations have changed. Concentrations change + * whenever the pressure or the mole fraction has changed. + * If it has changed, the recalculations should be done. + * + * Note this should be a lightweight function since it's + * part of all of the interfaces. + * + * @internal + */ + virtual bool update_C(); + + //! Update the temperature-dependent viscosity terms. + //! Updates the array of pure species viscosities, and the + //! weighting functions in the viscosity mixture rule. + /*! + * The flag m_visc_temp_ok is set to true. + */ + void updateViscosity_T(); + + //! Update the temperature-dependent parts of the mixture-averaged + //! thermal conductivity. + void updateCond_T(); + + //! Update the concentration parts of the viscosities + /*! + * Internal routine is run whenever the update_boolean + * is false. This routine will calculate + * internal values for the species viscosities. + * + * @internal + */ + void updateViscosities_C(); + + //! Update the binary diffusion coefficients wrt T. + /*! + * These are evaluated + * from the polynomial fits at unit pressure (1 Pa). + */ + void updateDiff_T(); + + + private: + + //! Number of species in the mixture + int m_nsp; + + //! Temperature dependence type + /*! + * The following coefficients are allowed to have simple + * temperature dependencies: + * mixture viscosity + * mixture thermal conductivity + * diffusitivy + * + * Types of temperature dependencies: + * 0 - Independent of temperature (only one implemented so far) + * 1 - extended arrhenius form + * 2 - power law form + */ + int tempDepType_; + + //! Composition dependence of the transport properties + /*! + * The following coefficients are allowed to have simple + * composition dependencies + * mixture viscosity + * mixture thermal conductivity + * + * + * Types of composition dependencies + * 0 - Solvent values (i.e., species 0) contributes only + * 1 - linear combination of mole fractions; + */ + int compositionDepType_; + + //! Minimum temperature applicable to the transport property eval + doublereal m_tmin; + + //! Maximum temperature applicable to the transport property evaluator + doublereal m_tmax; + + //! Local Copy of the molecular weights of the species + /*! + * Length is Equal to the number of species in the mechanism. + */ + vector_fp m_mw; + + //! Pure species viscosities in Arrhenius temperature-dependent form. + vector m_coeffVisc_Ns; + + //! Pure species thermal conductivities in Arrhenius temperature-dependent form. + /*! + * + */ + vector m_coeffLambda_Ns; + + + //! Pure species viscosities in Arrhenius temperature-dependent form. + vector m_coeffDiff_Ns; + + + //! Internal value of the gradient of the mole fraction vector + /*! + * Note, this is the only gradient value that can and perhaps + * should reflect the true state of the mole fractions in the + * application solution vector. In other words no cropping or + * massaging of the values to make sure they are above zero + * should occur. - developing .... + * + * m_nsp is the number of species in the fluid + * k is the species index + * n is the dimensional index (x, y, or z). It has a length + * equal to m_nDim + * + * m_Grad_X[n*m_nsp + k] + */ + vector_fp m_Grad_X; + + //! Internal value of the gradient of the Temperature vector + /*! + * Generally, if a transport property needs this + * in its evaluation it will look to this place + * to get it. + * + * No internal property is precalculated based on gradients. + * Gradients are assumed to be freshly updated before + * every property call. + */ + vector_fp m_Grad_T; + + //! Internal value of the gradient of the Pressure vector + /*! + * Generally, if a transport property needs this + * in its evaluation it will look to this place + * to get it. + * + * No internal property is precalculated based on gradients. + * Gradients are assumed to be freshly updated before + * every property call. + */ + vector_fp m_Grad_P; + + //! Internal value of the gradient of the Electric Voltage + /*! + * Generally, if a transport property needs this + * in its evaluation it will look to this place + * to get it. + * + * No internal property is precalculated based on gradients. + * Gradients are assumed to be freshly updated before + * every property call. + */ + vector_fp m_Grad_V; + + + // property values + + + + //! Vector of Species Diffusivities + /*! + * Depends on the temperature. We have set the pressure dependence + * to zero for this liquid phase constituitve model + * + * units m2/s + */ + vector_fp m_diffSpecies; + + //! Species viscosities + /*! + * Viscosity of the species + * Length = number of species + * + * Depends on the temperature. We have set the pressure dependence + * to zero for this model + * + * controlling update boolean -> m_visc_temp_ok + */ + vector_fp m_viscSpecies; + + //! Internal value of the species individual thermal conductivities + /*! + * Then a mixture rule is applied to get the solution conductivities + * + * Depends on the temperature and perhaps pressure, but + * not the species concentrations + * + * controlling update boolean -> m_cond_temp_ok + */ + vector_fp m_condSpecies; + + //! State of the mole fraction vector. + int m_iStateMF; + + //! Local copy of the mole fractions of the species in the phase + /*! + * The mole fractions here are assumed to be bounded by 0.0 and 1.0 + * and they are assumed to add up to one exactly. This mole + * fraction vector comes from the ThermoPhase object. Derivative + * quantities from this are referred to as bounded. + * + * Update info? + * length = m_nsp + */ + vector_fp m_molefracs; + + + //! Local copy of the concentrations of the species in the phase + /*! + * The concentrations are consistent with the m_molefracs + * vector which is bounded and sums to one. + * + * Update info? + * length = m_nsp + */ + vector_fp m_concentrations; + + //! Local copy of the total concentration. + /*! + * This is consistent with the m_concentrations[] and + * m_molefracs[] vector. + */ + doublereal concTot_; + + + + doublereal meanMolecularWeight_; + doublereal dens_; + + //! Local copy of the charge of each species + /*! + * Contains the charge of each species (length m_nsp) + */ + vector_fp m_chargeSpecies; + + + //! Current Temperature -> locally storred + /*! + * This is used to test whether new temperature computations + * should be performed. + */ + doublereal m_temp; + + + //! Current value of the pressure + doublereal m_press; + + + //! Saved value of the mixture thermal conductivity + doublereal m_lambda; + + //! Saved value of the mixture viscosity + doublereal m_viscmix; + + //! work space + /*! + * Length is equal to m_nsp + */ + vector_fp m_spwork; + + + + private: + //! Boolean indicating that the top-level mixture viscosity is current + /*! + * This is turned false for every change in T, P, or C. + */ + bool m_visc_mix_ok; + + //! Boolean indicating that weight factors wrt viscosity is current + bool m_visc_temp_ok; + + //! Boolean indicating that mixture diffusion coeffs are current + bool m_diff_mix_ok; + + //! Boolean indicating that binary diffusion coeffs are current + bool m_diff_temp_ok; + + //! Flag to indicate that the pure species conductivities + //! are current wrt the temperature + bool m_cond_temp_ok; + + //! Boolean indicating that mixture conductivity is current + bool m_cond_mix_ok; + + + //! Number of dimensions + /*! + * Either 1, 2, or 3 + */ + int m_nDim; + + private: + + //! Throw an exception if this method is invoked. + /*! + * This probably indicates something is not yet implemented. + * + * @pram msg Indicates the member function which is not implemented + */ + doublereal err(std::string msg) const; + + }; +} +#endif + + + + + + diff --git a/Cantera/src/transport/TransportBase.h b/Cantera/src/transport/TransportBase.h old mode 100755 new mode 100644 index dd0d7aae3..02b3d2c90 --- a/Cantera/src/transport/TransportBase.h +++ b/Cantera/src/transport/TransportBase.h @@ -45,6 +45,7 @@ namespace Cantera { const int cFtnTransport = 600; const int cLiquidTransport = 700; const int cAqueousTransport = 750; + const int cSimpleTransport = 770; const int cRadiativeTransport = 800; const int cWaterTransport = 721; @@ -159,6 +160,14 @@ namespace Cantera { virtual doublereal viscosity() { return err("viscosity"); } + //! Returns the pure species viscosities + /*! + * The units are Pa-s and the length is the number of species + * + * @param visc Vector of viscosities + */ + virtual void getSpeciesViscosities(doublereal* const visc) + { err("getSpeciesViscosities"); } /** * The bulk viscosity in Pa-s. The bulk viscosity is only @@ -169,9 +178,11 @@ namespace Cantera { virtual doublereal bulkViscosity() { return err("bulkViscosity"); } - - /** - * The thermal conductivity in W/m/K. + //! Returns the mixture thermal conductivity in W/m/K. + /*! + * Units are in W / m K or equivalently kg m / s3 K + * + * @return returns thermal conductivity in W/m/K. */ virtual doublereal thermalConductivity() { return err("thermalConductivity"); } @@ -273,9 +284,9 @@ namespace Cantera { * Get the mass fluxes [kg/m^2/s], given the thermodynamic * state at two nearby points. * @param state1 Array of temperature, density, and mass - * fractions for state 1. + * fractions for state 1. * @param state2 Array of temperature, density, and mass - * fractions for state 2. + * fractions for state 2. * @param delta Distance from state 1 to state 2 (m). */ virtual void getMassFluxes(const doublereal* state1, @@ -298,9 +309,13 @@ namespace Cantera { { err("getThermalDiffCoeffs"); } - /** - * Binary diffusion coefficients [m^2/s]. - */ + //! Returns the matrix of binary diffusion coefficients [m^2/s]. + /*! + * @param ld Inner stride for writing the two dimension diffusion + * coefficients into a one dimensional vector + * @param d Diffusion coefficient matrix (must be at least m_k * m_k + * in length. + */ virtual void getBinaryDiffCoeffs(const int ld, doublereal* const d) { err("getBinaryDiffCoeffs"); } @@ -348,20 +363,20 @@ namespace Cantera { /** * Called by TransportFactory to set parameters. */ - virtual bool init(TransportParams& tr) - { err("init"); return false; } + //virtual bool init(TransportParams& tr) + //{ err("init"); return false; } /** * Called by TransportFactory to set parameters. */ virtual bool initGas( GasTransportParams& tr ) - { err("init"); return false; } + { err("initGas"); return false; } /** * Called by TransportFactory to set parameters. */ virtual bool initLiquid( LiquidTransportParams& tr ) - { err("init"); return false; } + { err("initLiquid"); return false; } diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp old mode 100755 new mode 100644 index 2beffa8e9..bb3dfbcd9 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -20,6 +20,7 @@ #include "MixTransport.h" #include "SolidTransport.h" #include "DustyGasTransport.h" +#include "SimpleTransport.h" #ifdef WITH_IDEAL_SOLUTIONS #include "LiquidTransport.h" @@ -96,6 +97,21 @@ namespace Cantera { }; + /** + * getArrhenius() parses the xml element called Arrhenius. + * The Arrhenius expression is + * \f[ k = A T^(b) exp (-E_a / RT). \f] + */ + static void getArrhenius(const XML_Node& node, + doublereal& A, doublereal& b, doublereal& E) { + /* parse the children for the A, b, and E conponents. + */ + A = getFloat(node, "A", "toSI"); + b = getFloat(node, "b"); + E = getFloat(node, "E", "actEnergy"); + E /= GasConstant; + } + //////////////////// class TransportFactory methods ////////////// @@ -239,6 +255,7 @@ namespace Cantera { m_models["CK_Mix"] = CK_MixtureAveraged; m_models["Liquid"] = cLiquidTransport; m_models["Aqueous"] = cAqueousTransport; + m_models["Simple"] = cSimpleTransport; m_models["User"] = cUserTransport; m_models["None"] = None; //m_models["Radiative"] = cRadiative; @@ -319,6 +336,11 @@ namespace Cantera { dtr = (DustyGasTransport*)tr; dtr->initialize(phase, gastr); break; + case cSimpleTransport: + tr = new SimpleTransport(); + initLiquidTransport(tr, phase, log_level); + tr->setThermo(*phase); + break; #ifdef WITH_IDEAL_SOLUTIONS case cLiquidTransport: tr = new LiquidTransport; @@ -342,7 +364,7 @@ namespace Cantera { -/** + /** * Prepare to build a new kinetic-theory-based transport manager * for low-density gases. Uses polynomial fits to Monchick & Mason * collision integrals. @@ -473,13 +495,13 @@ namespace Cantera { -/** + /** * Prepare to build a new transport manager for liquids assuming that * viscosity transport data is provided in Arhennius form. */ void TransportFactory::setupLiquidTransport(std::ostream &flog, - const std::vector &transport_database, - thermo_t* thermo, int log_level, LiquidTransportParams& trParam) { + const std::vector &transport_database, + thermo_t* thermo, int log_level, LiquidTransportParams& trParam) { // constant mixture attributes trParam.thermo = thermo; @@ -506,7 +528,7 @@ namespace Cantera { XML_Node root, log; getLiquidTransportData(transport_database, log, - trParam.thermo->speciesNames(), trParam); + trParam.thermo->speciesNames(), trParam); //int i, j; //for (i = 0; i < nsp; i++) trParam.poly[i].resize(nsp); @@ -617,7 +639,7 @@ namespace Cantera { // set up Monchick and Mason collision integrals setupMM(flog, transport_database, thermo, mode, log_level, trParam); // do model-specific initialization - tran->init(trParam); + tran->initGas(trParam); #ifdef DEBUG_MODE if (m_verbose) { trParam.xml->XML_close(flog, "transport"); @@ -633,8 +655,8 @@ namespace Cantera { * class and calls setupLiquidTransport(). */ void TransportFactory::initLiquidTransport(Transport* tran, - thermo_t* thermo, - int log_level) { + thermo_t* thermo, + int log_level) { const std::vector & transport_database = thermo->speciesData(); @@ -651,7 +673,7 @@ namespace Cantera { #endif setupLiquidTransport(flog, transport_database, thermo, log_level, trParam); // do model-specific initialization - tran->init(trParam); + tran->initLiquid(trParam); #ifdef DEBUG_MODE if (m_verbose) { trParam.xml->XML_close(flog, "transport"); @@ -780,7 +802,7 @@ namespace Cantera { for (i = 0; i < nsp; i++) { const XML_Node& sp = *xspecies[i]; name = sp["name"]; - std::cout << "Processing node for " << name << std::endl; + // std::cout << "Processing node for " << name << std::endl; // put in a try block so that species with no 'transport' // child are skipped, instead of throwing an exception. @@ -883,13 +905,17 @@ namespace Cantera { const std::vector &names, LiquidTransportParams& trParam) { - string name; + std::string name; + /* + * Create a map of species names versus liquid transport data parameters + */ std::map datatable; doublereal A_visc, n_visc, Tact_visc, hydrodynamic_radius; doublereal A_thcond, n_thcond, Tact_thcond; + doublereal A_spdiff, n_spdiff, Tact_spdiff; int nsp = static_cast(xspecies.size()); - std::cout << "Size of xspecies " << nsp << std::endl; + std::cout << "Size of xspecies " << nsp << std::endl; // read all entries in database into 'datatable' and check for // errors. Note that this procedure validates all entries, not @@ -900,53 +926,189 @@ namespace Cantera { for (i = 0; i < nsp; i++) { const XML_Node& sp = *xspecies[i]; name = sp["name"]; - std::cout << "Processing node for " << name << std::endl; + vector_fp vCoeff; + // std::cout << "Processing node for " << name << std::endl; // put in a try block so that species with no 'transport' // child are skipped, instead of throwing an exception. try { - XML_Node& trNode = sp.child("transport"); + if (sp.hasChild("transport")) { + XML_Node& trNode = sp.child("transport"); - hydrodynamic_radius = getFloat(trNode, "hydrodynamic_radius"); + // Fill datatable with LiquidTransportData objects for error checking + // and then insertion into LiquidTransportData objects below. + LiquidTransportData data; + data.speciesName = name; - XML_Node& visc = trNode.child("viscosity"); - getArrhenius(visc, A_visc, n_visc, Tact_visc ); + /* + * hydrodynamic radius + * + * format: + * 3.0 + * 3.0 + */ + if (trNode.hasChild("hydrodynamic_radius")) { + XML_Node& hnode = trNode.child("hydrodynamic_radius"); + std::string model = lowercase(hnode["model"]); + if (model == "" || model == "constant") { + hydrodynamic_radius = hnode.fp_value(); + if (hydrodynamic_radius > 0.0) data.hydroradius = hydrodynamic_radius; + else throw TransportDBError(linenum, + "negative or zero hydrodynamic radius"); + data.model_hydroradius = LTR_MODEL_CONSTANT; + } else { + throw CanteraError(" TransportFactory::getLiquidTransportData", + "Unknown model for hydrodynamic_radius:" + model); + } + } - XML_Node& thermCond = trNode.child("thermal_conductivity"); - getArrhenius(thermCond, A_thcond, n_thcond, Tact_thcond ); + /* + * viscosity + * + * format: + * 3.0 + * 3.0 + * + * 1.0 + * 2.0 + * 3.0 + * + * + * + * 0.0. 1.0, 2.0, 3.0, 4.0 + * + * + */ + if (trNode.hasChild("viscosity")) { + XML_Node& vnode = trNode.child("viscosity"); + std::string model = lowercase(vnode["model"]); + if (model == "" || model == "constant") { + A_visc = vnode.fp_value(); + if (A_visc > 0.0) (data.viscCoeffs).push_back(A_visc); + else throw TransportDBError(linenum, + "negative or zero viscosity"); + data.model_viscosity = LTR_MODEL_CONSTANT; + } else if (model == "arrhenius") { + getArrhenius(vnode, A_visc, n_visc, Tact_visc); + if (A_visc <= 0.0) { + throw TransportDBError(linenum, "negative or zero viscosity"); + } + (data.viscCoeffs).push_back(A_visc); + (data.viscCoeffs).push_back(n_visc); + (data.viscCoeffs).push_back(Tact_visc); + data.model_viscosity = LTR_MODEL_ARRHENIUS; + } else if (model == "coeff") { + getFloatArray(vnode, vCoeff, true); + data.viscCoeffs = vCoeff; + vCoeff.clear(); + data.model_viscosity = LTR_MODEL_COEFF; + } else { + throw CanteraError(" TransportFactory::getLiquidTransportData", + "Unknown model for viscosity:" + vnode["model"]); + } + } - // Fill datatable with LiquidTransportData objects for error checking - // and then insertion into LiquidTransportData objects below. - LiquidTransportData data; - data.speciesName = name; + /* + * thermal_conductivity + * + * format: + * 3.0 + * 3.0 + * + * 1.0 + * 2.0 + * 3.0 + * + * + * + * 0.0. 1.0, 2.0, 3.0, 4.0 + * + * + */ + if (trNode.hasChild("thermal_conductivity")) { + XML_Node& tnode = trNode.child("thermal_conductivity"); + std::string model = lowercase(tnode["model"]); + if (model == "" || model == "constant") { + A_thcond = tnode.fp_value(); + if (A_thcond > 0.0) (data.thermalCondCoeffs).push_back(A_thcond); + else throw TransportDBError(linenum, + "negative or zero thermal_conductivity"); + data.model_thermalCond = LTR_MODEL_CONSTANT; + } else if (model == "arrhenius") { + getArrhenius(tnode, A_thcond, n_thcond, Tact_thcond); + if (A_thcond <= 0.0) { + throw TransportDBError(linenum, "negative or zero thermal_conductivity"); + } + (data.thermalCondCoeffs).push_back(A_thcond); + (data.thermalCondCoeffs).push_back(n_thcond); + (data.thermalCondCoeffs).push_back(Tact_thcond); + data.model_thermalCond = LTR_MODEL_ARRHENIUS; + } else if (model == "coeff") { + getFloatArray(tnode, vCoeff, true); + data.thermalCondCoeffs = vCoeff; + vCoeff.clear(); + data.model_thermalCond = LTR_MODEL_COEFF; + } else { + throw CanteraError(" TransportFactory::getLiquidTransportData", + "Unknown model for thermal_conductivity:" + tnode["model"]); + } + } - if ( hydrodynamic_radius > 0.0) data.hydroradius = hydrodynamic_radius; - else throw TransportDBError(linenum, - "negative or zero hydrodynamic radius"); + + /* + * speciesDiffusivity + * + * format: + * 3.0 + * 3.0 + * + * 1.0 + * 2.0 + * 3.0 + * + * + * + * 0.0. 1.0, 2.0, 3.0, 4.0 + * + * + */ + if (trNode.hasChild("speciesDiffusivity")) { + XML_Node& dnode = trNode.child("speciesDiffusivity"); + std::string model = lowercase(dnode["model"]); + if (model == "" || model == "constant") { + A_spdiff = dnode.fp_value(); + if (A_spdiff > 0.0) (data.speciesDiffusivityCoeffs).push_back(A_spdiff); + else throw TransportDBError(linenum, + "negative or zero speciesDiffusivity"); + data.model_speciesDiffusivity = LTR_MODEL_CONSTANT; + } else if (model == "arrhenius") { + getArrhenius(dnode, A_spdiff, n_spdiff, Tact_spdiff); + if (A_spdiff <= 0.0) { + throw TransportDBError(linenum, "negative or zero speciesDiffusivity"); + } + (data.speciesDiffusivityCoeffs).push_back(A_spdiff); + (data.speciesDiffusivityCoeffs).push_back(n_spdiff); + (data.speciesDiffusivityCoeffs).push_back(Tact_spdiff); + data.model_speciesDiffusivity = LTR_MODEL_ARRHENIUS; + } else if (model == "coeff") { + getFloatArray(dnode, vCoeff, true); + data.speciesDiffusivityCoeffs = vCoeff; + data.model_speciesDiffusivity = LTR_MODEL_COEFF; + } else { + throw CanteraError(" TransportFactory::getLiquidTransportData", + "Unknown model for speciesDiffusivity:" + dnode["model"]); + } + } - if (A_visc >= 0.0) { - data.viscCoeffs[0] = A_visc; - data.viscCoeffs[1] = n_visc; - data.viscCoeffs[2] = Tact_visc; + datatable[name] = data; } - else throw TransportDBError(linenum, - "negative pre-exponential for viscosity"); - - if (A_thcond >= 0.0) { - data.thermalCondCoeffs[0] = A_thcond; - data.thermalCondCoeffs[1] = n_thcond; - data.thermalCondCoeffs[2] = Tact_thcond; - } - else throw TransportDBError(linenum, - "negative pre-exponential for thermalCondoctivity"); - - datatable[name] = data; } catch(CanteraError) { ; } } + trParam.LTData.clear(); for (i = 0; i < trParam.nsp_; i++) { LiquidTransportData& trdat = datatable[names[i]]; @@ -960,24 +1122,36 @@ namespace Cantera { } // parameters should be converted to SI units before storing - - trParam.visc_A[i] = trdat.viscCoeffs[0] ; - trParam.visc_n[i] = trdat.viscCoeffs[1] ; - trParam.visc_Tact[i] = trdat.viscCoeffs[2] ; + if (trdat.viscCoeffs.size() > 0) { + trParam.visc_A[i] = trdat.viscCoeffs[0] ; + } + if (trdat.viscCoeffs.size() > 2) { + trParam.visc_n[i] = trdat.viscCoeffs[1] ; + trParam.visc_Tact[i] = trdat.viscCoeffs[2] ; + } - trParam.thermCond_A[i] = trdat.thermalCondCoeffs[0] ; - trParam.thermCond_n[i] = trdat.thermalCondCoeffs[1] ; - trParam.thermCond_Tact[i] = trdat.thermalCondCoeffs[2] ; + if (trdat.thermalCondCoeffs.size() > 0) { + trParam.thermCond_A[i] = trdat.thermalCondCoeffs[0] ; + } + if (trdat.thermalCondCoeffs.size() > 2) { + trParam.thermCond_n[i] = trdat.thermalCondCoeffs[1] ; + trParam.thermCond_Tact[i] = trdat.thermalCondCoeffs[2] ; + } // Angstroms -> meters trParam.hydroRadius[i] = 1.e-10 * trdat.hydroradius; + /* + * this is a much more general way to handle the transfer + * -> calling the default copy constructor for LiquidTransportData + */ + trParam.LTData.push_back(trdat); } - // Need to identify a method to obtain interaction matrices. - // This will fill LiquidTransportParams members visc_Eij, visc_Sij - trParam.visc_Eij.resize(trParam.nsp_,trParam.nsp_); - cout << "No support for species viscosity interactions in TransportFactory.cpp" << endl; + // Need to identify a method to obtain interaction matrices. + // This will fill LiquidTransportParams members visc_Eij, visc_Sij + trParam.visc_Eij.resize(trParam.nsp_,trParam.nsp_); + cout << "No support for species viscosity interactions in TransportFactory.cpp" << endl; } diff --git a/Cantera/src/transport/TransportFactory.h b/Cantera/src/transport/TransportFactory.h index 489ef7156..c89f1b7d1 100755 --- a/Cantera/src/transport/TransportFactory.h +++ b/Cantera/src/transport/TransportFactory.h @@ -34,6 +34,7 @@ #include "ct_defs.h" #include "TransportBase.h" #include "FactoryBase.h" +#include "LiquidTransportData.h" #if defined(THREAD_SAFE_CANTERA) #include @@ -61,20 +62,12 @@ namespace Cantera { doublereal rotRelaxNumber; }; - struct LiquidTransportData { - LiquidTransportData() : speciesName("-"), - hydroradius(-1) {} - std::string speciesName; - doublereal hydroradius; - vector_fp viscCoeffs; - vector_fp thermalCondCoeffs; - }; // forward references class MMCollisionInt; class GasTransportParams; class LiquidTransportParams; - class XML_Node; + class XML_Node; //! The purpose of TransportFactory is to create new instances of @@ -173,9 +166,19 @@ namespace Cantera { XML_Node& log, const std::vector& names, GasTransportParams& tr); + + //! Read transport property data from a file for a list of species. + /*! + * + * Given the name of a file containing transport property + * parameters and a list of species names, this method returns an + * instance of TransportParams containing the transport data for + * these species read from the file. + * + */ void getLiquidTransportData(const std::vector &db, - XML_Node& log, const std::vector& names, - LiquidTransportParams& tr); + XML_Node& log, const std::vector& names, + LiquidTransportParams& tr); /** Generate polynomial fits to viscosity, conductivity, and * binary diffusion coefficients */ @@ -193,8 +196,8 @@ namespace Cantera { void setupLiquidTransport(std::ostream &flog, const std::vector &transport_database, - thermo_t* thermo, int log_level, - LiquidTransportParams& tr); + thermo_t* thermo, int log_level, + LiquidTransportParams& tr); /// Second-order correction to the binary diffusion coefficients @@ -208,21 +211,6 @@ namespace Cantera { const GasTransportParams& tr, doublereal& f_eps, doublereal& f_sigma); - /** - * getArrhenius() parses the xml element called Arrhenius. - * The Arrhenius expression is - * \f[ k = A T^(b) exp (-E_a / RT). \f] - */ - static void getArrhenius(const XML_Node& node, - doublereal& A, doublereal& b, doublereal& E) { - /* parse the children for the A, b, and E conponents. - */ - A = getFloat(node, "A", "toSI"); - b = getFloat(node, "b"); - E = getFloat(node, "E", "actEnergy"); - E /= GasConstant; - } - //! Boolean indicating whether to turn on verbose printing bool m_verbose;