From 0127285b881ad2e80a00206a5bfb4a1b2fa7ce0a Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 14 Feb 2009 19:33:52 +0000 Subject: [PATCH] Added stubs for liquid phase and aqueous phase transport. --- Cantera/src/transport/AqueousTransport.cpp | 735 ++++++++++++++++++ Cantera/src/transport/AqueousTransport.h | 621 +++++++++++++++ Cantera/src/transport/LiquidTransport.cpp | 722 +++++++++++++++++ Cantera/src/transport/LiquidTransport.h | 634 +++++++++++++++ Cantera/src/transport/LiquidTransportParams.h | 62 ++ Cantera/src/transport/Makefile.in | 5 +- Cantera/src/transport/MixTransport.h | 8 + Cantera/src/transport/TransportBase.h | 1 + Cantera/src/transport/TransportFactory.cpp | 14 +- 9 files changed, 2799 insertions(+), 3 deletions(-) create mode 100644 Cantera/src/transport/AqueousTransport.cpp create mode 100644 Cantera/src/transport/AqueousTransport.h create mode 100644 Cantera/src/transport/LiquidTransport.cpp create mode 100644 Cantera/src/transport/LiquidTransport.h create mode 100644 Cantera/src/transport/LiquidTransportParams.h diff --git a/Cantera/src/transport/AqueousTransport.cpp b/Cantera/src/transport/AqueousTransport.cpp new file mode 100644 index 000000000..085b07391 --- /dev/null +++ b/Cantera/src/transport/AqueousTransport.cpp @@ -0,0 +1,735 @@ +/** + * @file MixTransport.cpp + * Mixture-averaged transport properties for ideal gas mixtures. + */ +/* + * $Revision$ + * $Date$ + */ + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ThermoPhase.h" +#include "LiquidTransport.h" + +#include "utilities.h" +#include "TransportParams.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-20 + + +namespace Cantera { + + //////////////////// class AqueousTransport methods ////////////// + + + AqueousTransport::AqueousTransport() : + m_nsp(0), + m_tmin(-1.0), + m_tmax(100000.), + m_iStateMF(-1), + m_temp(-1.0), + m_logt(0.0), + m_sqrt_t(-1.0), + m_t14(-1.0), + m_t32(-1.0), + m_sqrt_kbt(-1.0), + m_press(-1.0), + m_lambda(-1.0), + m_viscmix(-1.0), + m_viscmix_ok(false), + m_viscwt_ok(false), + m_spvisc_ok(false), + m_diffmix_ok(false), + m_bindiff_ok(false), + m_spcond_ok(false), + m_condmix_ok(false), + m_mode(-1000), + m_debug(false), + m_nDim(1) + { + + + } + + + // Initialize the object + /* + * This is where we dimension everything. + */ + bool AqueousTransport::init(TransportParams& 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()); + + // copy polynomials and parameters into local storage + m_poly = tr.poly; + m_visccoeffs = tr.visccoeffs; + m_condcoeffs = tr.condcoeffs; + m_diffcoeffs = tr.diffcoeffs; + + m_mode = tr.mode; + m_diam = tr.diam; + m_eps = tr.eps; + m_alpha = tr.alpha; + + m_phi.resize(m_nsp, m_nsp, 0.0); + + + m_wratjk.resize(m_nsp, m_nsp, 0.0); + m_wratkj1.resize(m_nsp, m_nsp, 0.0); + int j, k; + for (j = 0; j < m_nsp; j++) + for (k = j; k < m_nsp; k++) { + m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]); + m_wratjk(k,j) = sqrt(m_wratjk(j,k)); + m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]); + } + + m_polytempvec.resize(5); + m_visc.resize(m_nsp); + m_sqvisc.resize(m_nsp); + m_cond.resize(m_nsp); + m_bdiff.resize(m_nsp, 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); + m_Grad_mu.resize(m_nDim * m_nsp, 0.0); + + + // set all flags to false + m_viscmix_ok = false; + m_viscwt_ok = false; + m_spvisc_ok = false; + m_spcond_ok = false; + m_condmix_ok = false; + m_spcond_ok = false; + m_diffmix_ok = false; + + return true; + } + + + /********************************************************* + * + * Public methods + * + *********************************************************/ + + + /****************** viscosity ******************************/ + + /* + * The viscosity is computed using the Wilke mixture rule. + * \f[ + * \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}. + * \f] + * Here \f$ \mu_k \f$ is the viscosity of pure species \e k, + * and + * \f[ + * \Phi_{k,j} = \frac{\left[1 + * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2} + * {\sqrt{8}\sqrt{1 + M_k/M_j}} + * \f] + * @see updateViscosity_T(); + */ + doublereal AqueousTransport::viscosity() { + + update_T(); + update_C(); + + if (m_viscmix_ok) return m_viscmix; + + // update m_visc[] and m_phi[] if necessary + if (!m_viscwt_ok) updateViscosity_T(); + + multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork)); + + m_viscmix = 0.0; + for (int k = 0; k < m_nsp; k++) { + m_viscmix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom; + } + return m_viscmix; + } + + + /******************* binary diffusion coefficients **************/ + + + void AqueousTransport::getBinaryDiffCoeffs(int ld, doublereal* d) { + int i,j; + + update_T(); + + // if necessary, evaluate the binary diffusion coefficents + // from the polynomial fits + if (!m_bindiff_ok) updateDiff_T(); + doublereal pres = m_thermo->pressure(); + + doublereal rp = 1.0/pres; + for (i = 0; i < m_nsp; i++) + for (j = 0; j < m_nsp; j++) { + d[ld*j + i] = rp * m_bdiff(i,j); + } + } + + + void AqueousTransport::getMobilities(doublereal* 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 AqueousTransport::set_Grad_V(const doublereal* grad_V) { + for (int a = 0; a < m_nDim; a++) { + m_Grad_V[a] = grad_V[a]; + } + } + + void AqueousTransport::set_Grad_T(const doublereal* grad_T) { + for (int a = 0; a < m_nDim; a++) { + m_Grad_T[a] = grad_T[a]; + } + } + + void AqueousTransport::set_Grad_X(const doublereal* grad_X) { + int itop = m_nDim * m_nsp; + for (int i = 0; i < itop; i++) { + m_Grad_X[i] = grad_X[i]; + } + } + + + /****************** thermal conductivity **********************/ + + /* + * The thermal conductivity is computed from the following mixture rule: + * \[ + * \lambda = 0.5 \left( \sum_k X_k \lambda_k + * + \frac{1}{\sum_k X_k/\lambda_k}\right) + * \] + */ + doublereal AqueousTransport::thermalConductivity() { + int k; + + update_T(); + update_C(); + + if (!m_spcond_ok) updateCond_T(); + if (!m_condmix_ok) { + doublereal sum1 = 0.0, sum2 = 0.0; + for (k = 0; k < m_nsp; k++) { + sum1 += m_molefracs[k] * m_cond[k]; + sum2 += m_molefracs[k] / m_cond[k]; + } + m_lambda = 0.5*(sum1 + 1.0/sum2); + } + return m_lambda; + } + + + /****************** thermal diffusion coefficients ************/ + + /** + * 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 AqueousTransport::getThermalDiffCoeffs(doublereal* dt) { + int k; + for (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 AqueousTransport::getSpeciesFluxes(int ndim, + 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); + } + + /** + * @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 AqueousTransport::getSpeciesFluxesExt(int ldf, doublereal* fluxes) { + int n, k; + + 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]. + * + * For the single species case or the pure fluid case + * the routine returns the self-diffusion coefficient. + * This is need to avoid a Nan result in the formula + * below. + */ + void AqueousTransport::getMixDiffCoeffs(doublereal* d) { + + update_T(); + update_C(); + + // update the binary diffusion coefficients if necessary + if (!m_bindiff_ok) updateDiff_T(); + + int k, j; + doublereal mmw = m_thermo->meanMolecularWeight(); + doublereal sumxw = 0.0, sum2; + doublereal p = m_press; + if (m_nsp == 1) { + d[0] = m_bdiff(0,0) / p; + } else { + for (k = 0; k < m_nsp; k++) sumxw += m_molefracs[k] * m_mw[k]; + for (k = 0; k < m_nsp; k++) { + sum2 = 0.0; + for (j = 0; j < m_nsp; j++) { + if (j != k) { + sum2 += m_molefracs[j] / m_bdiff(j,k); + } + } + if (sum2 <= 0.0) { + d[k] = m_bdiff(k,k) / p; + } else { + d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2); + } + } + } + } + + + // 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 + */ + void AqueousTransport::update_T() + { + doublereal t = m_thermo->temperature(); + if (t == m_temp) return; + if (t < 0.0) { + throw CanteraError("AqueousTransport::update_T", + "negative temperature "+fp2str(t)); + } + + // Compute various functions of temperature + m_temp = t; + m_logt = log(m_temp); + m_kbt = Boltzmann * m_temp; + m_sqrt_t = sqrt(m_temp); + m_t14 = sqrt(m_sqrt_t); + m_t32 = m_temp * m_sqrt_t; + m_sqrt_kbt = sqrt(Boltzmann*m_temp); + + // compute powers of log(T) + m_polytempvec[0] = 1.0; + m_polytempvec[1] = m_logt; + m_polytempvec[2] = m_logt*m_logt; + m_polytempvec[3] = m_logt*m_logt*m_logt; + m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt; + + // temperature has changed, so polynomial temperature + // interpolations will need to be reevaluated. + // Set all of these flags to false + m_viscmix_ok = false; + m_spvisc_ok = false; + m_viscwt_ok = false; + m_spcond_ok = false; + m_diffmix_ok = false; + m_bindiff_ok = false; + m_condmix_ok = false; + + // For now, for a concentration redo also + m_iStateMF = -1; + } + + /** + * @internal This is called the first time any transport property + * is requested from Mixture after the concentrations + * have changed. + */ + void AqueousTransport::update_C() + { + + doublereal pres = m_thermo->pressure(); + // Check for changes in the mole fraction vector. + int iStateNew = m_thermo->getIStateMF(); + if (iStateNew == m_iStateMF) { + if (pres == m_press) { + return; + } + } else { + m_iStateMF = iStateNew; + } + m_press = pres; + + // signal that concentration-dependent quantities will need to + // be recomputed before use, and update the local mole + // fractions. + + m_viscmix_ok = false; + m_diffmix_ok = false; + m_condmix_ok = false; + + m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); + + // add an offset to avoid a pure species condition or + // negative mole fractions. MIN_X is 1.0E-20, a value + // which is below the additive machine precision of mole fractions. + int k; + for (k = 0; k < m_nsp; k++) { + m_molefracs[k] = fmaxx(MIN_X, m_molefracs[k]); + } + } + + + /************************************************************************* + * + * methods to update temperature-dependent properties + * + *************************************************************************/ + + /** + * Update the temperature-dependent parts of the mixture-averaged + * thermal conductivity. + */ + void AqueousTransport::updateCond_T() { + + int k; + if (m_mode == CK_Mode) { + for (k = 0; k < m_nsp; k++) { + m_cond[k] = exp(dot4(m_polytempvec, m_condcoeffs[k])); + } + } + else { + for (k = 0; k < m_nsp; k++) { + m_cond[k] = m_sqrt_t*dot5(m_polytempvec, m_condcoeffs[k]); + } + } + m_spcond_ok = true; + m_condmix_ok = false; + } + + + /** + * Update the binary diffusion coefficients. These are evaluated + * from the polynomial fits at unit pressure (1 Pa). + */ + void AqueousTransport::updateDiff_T() { + + // evaluate binary diffusion coefficients at unit pressure + int i,j; + int ic = 0; + if (m_mode == CK_Mode) { + for (i = 0; i < m_nsp; i++) { + for (j = i; j < m_nsp; j++) { + m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic])); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } + else { + for (i = 0; i < m_nsp; i++) { + for (j = i; j < m_nsp; j++) { + m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec, + m_diffcoeffs[ic]); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } + + m_bindiff_ok = true; + m_diffmix_ok = false; + } + + + /** + * Update the pure-species viscosities. + */ + void AqueousTransport::updateSpeciesViscosities() { + + int k; + if (m_mode == CK_Mode) { + for (k = 0; k < m_nsp; k++) { + m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k])); + m_sqvisc[k] = sqrt(m_visc[k]); + } + } + else { + for (k = 0; k < m_nsp; k++) { + // the polynomial fit is done for sqrt(visc/sqrt(T)) + m_sqvisc[k] = m_t14*dot5(m_polytempvec, m_visccoeffs[k]); + m_visc[k] = (m_sqvisc[k]*m_sqvisc[k]); + } + } + m_spvisc_ok = true; + } + + + /** + * 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 AqueousTransport::updateViscosity_T() { + doublereal vratiokj, wratiojk, factor1; + + if (!m_spvisc_ok) updateSpeciesViscosities(); + + // see Eq. (9-5.15) of Reid, Prausnitz, and Poling + int j, k; + for (j = 0; j < m_nsp; j++) { + for (k = j; k < m_nsp; k++) { + vratiokj = m_visc[k]/m_visc[j]; + wratiojk = m_mw[j]/m_mw[k]; + + // Note that m_wratjk(k,j) holds the square root of + // m_wratjk(j,k)! + factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j); + m_phi(k,j) = factor1*factor1 / + (SqrtEight * m_wratkj1(j,k)); + m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk); + } + } + m_viscwt_ok = true; + } + + /** + * This function returns a Transport data object for a given species. + * + */ + struct GasTransportData AqueousTransport:: + getGasTransportData(int kSpecies) + { + struct GasTransportData td; + td.speciesName = m_thermo->speciesName(kSpecies); + + + td.wellDepth = m_eps[kSpecies] / Boltzmann; + td.diameter = m_diam(kSpecies, kSpecies) * 1.0E10; + td.polarizability = m_alpha[kSpecies] * 1.0E30; + + + return td; + } + + + /* + * + * Solve for the diffusional velocities in the Stefan-Maxwell equations + * + */ + void AqueousTransport::stefan_maxwell_solve() { + int i, j, a; + + int VIM = 2; + m_B.resize(m_nsp, VIM); + //! grab a local copy of the molecular weights + const vector_fp& M = m_thermo->molecularWeights(); + + + //! get the mean molecular weight of the mixture + //double M_mix = m_thermo->meanMolecularWeight(); + + + //! get the concentration of the mixture + //double rho = m_thermo->density(); + //double c = rho/M_mix; + + + m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); + + double T = m_thermo->temperature(); + + + /* electrochemical potential gradient */ + for (i = 0; i < m_nsp; i++) { + for (a = 0; a < VIM; a++) { + m_Grad_mu[a*m_nsp + i] = m_chargeSpecies[i] * Faraday * m_Grad_V[a] + + (GasConstant*T/m_molefracs[i]) * m_Grad_X[a*m_nsp+i]; + } + } + + /* + * Just for Note, m_A(i,j) refers to the ith row and jth column. + * They are still fortran ordered, so that i varies fastest. + */ + switch ( VIM ) { + case 1: /* 1-D approximation */ + m_B(0,0) = 0.0; + for (j = 0; j < m_nsp; j++) { + m_A(0,j) = 1.0; + } + for (i = 1; i < m_nsp; i++){ + m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T); + for (j = 0; j < m_nsp; j++){ + if (j != i) { + m_A(i,j) = m_molefracs[i] / ( M[j] * m_DiffCoeff_StefMax(i,j)); + m_A(i,i) -= m_molefracs[j] / ( M[i] * m_DiffCoeff_StefMax(i,j)); + } + else if (j == i) { + m_A(i,i) = 0.0; + } + } + } + + //! invert and solve the system Ax = b. Answer is in m_B + solve(m_A, m_B.ptrColumn(0)); + + m_flux = m_B; + + + break; + case 2: /* 2-D approximation */ + m_B(0,0) = 0.0; + m_B(0,1) = 0.0; + for (j = 0; j < m_nsp; j++) { + m_A(0,j) = 1.0; + } + for (i = 1; i < m_nsp; i++){ + m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T); + m_B(i,1) = m_concentrations[i] * m_Grad_mu[m_nsp + i] / (GasConstant * T); + for (j = 0; j < m_nsp; j++){ + if (j != i) { + m_A(i,j) = m_molefracs[i] / ( M[j] * m_DiffCoeff_StefMax(i,j)); + m_A(i,i) -= m_molefracs[j] / ( M[i] * m_DiffCoeff_StefMax(i,j)); + } + else if (j == i) { + m_A(i,i) = 0.0; + } + } + } + + //! invert and solve the system Ax = b. Answer is in m_B + solve(m_A, m_B); + + m_flux = m_B; + + + break; + + case 3: /* 3-D approximation */ + m_B(0,0) = 0.0; + m_B(0,1) = 0.0; + m_B(0,2) = 0.0; + for (j = 0; j < m_nsp; j++) { + m_A(0,j) = 1.0; + } + for (i = 1; i < m_nsp; i++){ + m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T); + m_B(i,1) = m_concentrations[i] * m_Grad_mu[m_nsp + i] / (GasConstant * T); + m_B(i,2) = m_concentrations[i] * m_Grad_mu[2*m_nsp + i] / (GasConstant * T); + for (j = 0; j < m_nsp; j++){ + if (j != i) { + m_A(i,j) = m_molefracs[i] / ( M[j] * m_DiffCoeff_StefMax(i,j)); + m_A(i,i) -= m_molefracs[j] / ( M[i] * m_DiffCoeff_StefMax(i,j)); + } + else if (j == i) { + m_A(i,i) = 0.0; + } + } + } + + //! invert and solve the system Ax = b. Answer is in m_B + solve(m_A, m_B); + + m_flux = m_B; + + + break; + default: + printf("uninmplemetnd\n"); + throw CanteraError("routine", "not done"); + break; + } + + + } +} diff --git a/Cantera/src/transport/AqueousTransport.h b/Cantera/src/transport/AqueousTransport.h new file mode 100644 index 000000000..3b2e68dbb --- /dev/null +++ b/Cantera/src/transport/AqueousTransport.h @@ -0,0 +1,621 @@ +/** + * @file LiquidTransport.h + * Header file defining class AqueousTransport + */ +/* + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_AQUEOUSTRAN_H +#define CT_AQYEOUSTRAN_H + +using namespace std; + +// Cantera includes +#include "TransportBase.h" +#include "DenseMatrix.h" + + +#include +#include +#include +#include +#include + +namespace Cantera { + + + class TransportParams; + + + //! Class AqueousTransport 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 AqueousTransport : public Transport { + + public: + + //! virtual destructor + virtual ~AqueousTransport() {} + + //! Return the model id for this transport parameterization + virtual int model() { return cAqueousTransport; } + + //! overloaded base class methods + + //! Returns the viscosity of the solution + /*! + * The viscosity is computed using the Wilke mixture rule. + * \f[ + * \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}. + * \f] + * Here \f$ \mu_k \f$ is the viscosity of pure species \e k, + * and + * \f[ + * \Phi_{k,j} = \frac{\left[1 + * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2} + * {\sqrt{8}\sqrt{1 + M_k/M_j}} + * \f] + * @see updateViscosity_T(); + * + * Controlling update boolean m_viscmix_ok + */ + virtual doublereal viscosity(); + + //! Returns the pure species viscosities + /*! + * + * Controlling update boolean = m_viscwt_ok + */ + virtual void getSpeciesViscosities(doublereal* visc) + { updateViscosity_T(); copy(m_visc.begin(), m_visc.end(), visc); } + + virtual void getThermalDiffCoeffs(doublereal* dt); + + //! Return the thermal conductivity of the solution + /*! + * The thermal conductivity is computed from the following mixture rule: + * \f[ + * \lambda = 0.5 \left( \sum_k X_k \lambda_k + * + \frac{1}{\sum_k X_k/\lambda_k}\right) + * \f] + * + * Controlling update boolean = m_condmix_ok + */ + virtual doublereal thermalConductivity(); + + //! Returns the binary diffusion coefficients + /*! + * @param ld + * @param d + */ + virtual void getBinaryDiffCoeffs(int ld, doublereal* 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* d); + + + //! Get the Mobilities + /*! + * @param mobil + */ + virtual void getMobilities(doublereal* mobil); + + //! Specify the value of the gradient of the voltage + /*! + * + * @param grad_V Gradient of the voltage (length num dimensions); + */ + virtual void set_Grad_V(const doublereal* 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* 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* grad_X); + + //! 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 + */ + virtual void update_T(); + + //! Handles the effects of changes in the mixture concentration + /*! + * This is called the first time any transport property + * is requested from Mixture after the concentrations + * have changed. + * + * @internal + */ + virtual void update_C(); + + /** + * @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, + doublereal* grad_T, + int ldx, const doublereal* grad_X, + int ldf, doublereal* 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 getSpeciesFluxesExt(int ldf, doublereal* fluxes); + + + //! 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 init(TransportParams& tr); + + friend class TransportFactory; + + /** + * Return a structure containing all of the pertinent parameters + * about a species that was used to construct the Transport + * properties in this object. + * + * @param k Species number to obtain the properties about. + */ + struct GasTransportData getGasTransportData(int k); + + + //! Solve the stefan_maxell equations for the diffusive fluxes. + void stefan_maxwell_solve(); + protected: + + //! default constructor + AqueousTransport(); + + private: + + + + //! Number of species in the mixture + int m_nsp; + + //! 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; + + // polynomial fits + vector > m_poly; + + //! Polynomial coefficients of the viscosity + /*! + * These express the temperature dependendence of the pures + * species viscosities. + */ + vector m_visccoeffs; + + //! Polynomial coefficients of the conductivities + /*! + * These express the temperature dependendence of the pures + * species conductivities + */ + vector m_condcoeffs; + + //! Polynomial coefficients of the binary diffusion coefficients + /*! + * These express the temperature dependendence of the + * binary diffusivities. An overall pressure dependence is then + * added. + */ + vector m_diffcoeffs; + + + //! Internal value of the gradient of the mole fraction vector + /*! + * 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 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; + + //! Gradient of the electrochemical potential + /*! + * m_nsp is the number of species in the fluid + * k is the species index + * n is the dimensional index (x, y, or z) + * + * m_Grad_mu[n*m_nsp + k] + */ + vector_fp m_Grad_mu; + + // property values + + //! Array of Binary Diffusivities + /*! + * This has a size equal to nsp x nsp + * It is a symmetric matrix. + * D_ii is undefined. + * + * units m2/sec + */ + DenseMatrix m_bdiff; + + //! Species viscosities + /*! + * Viscosity of the species + * Length = number of species + * + * Depends on the temperature and perhaps pressure, but + * not the species concentrations + * + * controlling update boolean -> m_spvisc_ok + */ + vector_fp m_visc; + + //! Sqrt of the species viscosities + /*! + * The sqrt(visc) is used in the mixing formulas + * Length = m_nsp + * + * Depends on the temperature and perhaps pressure, but + * not the species concentrations + * + * controlling update boolean m_spvisc_ok + */ + vector_fp m_sqvisc; + + //! 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_spcond_ok + */ + vector_fp m_cond; + + //! Polynomials of the log of the temperature + vector_fp m_polytempvec; + + //! State of the mole fraction vector. + int m_iStateMF; + + //! Local copy of the mole fractions of the species in the phase + /*! + * Update info? + * length = m_nsp + */ + vector_fp m_molefracs; + + //! Local copy of the concentrations of the species in the phase + /*! + * Update info? + * length = m_nsp + */ + vector_fp m_concentrations; + + //! Local copy of the charge of each species + /*! + * Contains the charge of each species (length m_nsp) + */ + vector_fp m_chargeSpecies; + + //! Stefan-Maxwell Diffusion Coefficients at T, P and C + /*! + * These diffusion coefficients are considered to be + * a function of Temperature, Pressure, and Concentration. + */ + DenseMatrix m_DiffCoeff_StefMax; + + //! viscosity weighting functions + DenseMatrix m_phi; + + //! Matrix of the ratios of the species molecular weights + /*! + * m_wratjk(i,j) = (m_mw[j]/m_mw[k])**0.25 + */ + DenseMatrix m_wratjk; + + //! Matrix of the ratios of the species molecular weights + /*! + * m_wratkj1(i,j) = (1.0 + m_mw[k]/m_mw[j])**0.5 + */ + DenseMatrix m_wratkj1; + + //! RHS to the stefan-maxwell equation + Array2D m_B; + + //! Matrix for the stefan maxwell equation. + DenseMatrix m_A; + + //! Internal storage for the species LJ well depth + vector_fp m_eps; + + //! Internal storage for species polarizability + vector_fp m_alpha; + + //! Current Temperature -> locally storred + /*! + * This is used to test whether new temperature computations + * should be performed. + */ + doublereal m_temp; + + //! Current log(T) + doublereal m_logt; + + //! Current value of kT + doublereal m_kbt; + + //! Current Temperature **0.5 + doublereal m_sqrt_t; + + //! Current Temperature **0.25 + doublereal m_t14; + + //! Current Temperature **1.5 + doublereal m_t32; + + //! Current temperature function + /*! + * This is equal to sqrt(Boltzmann * T) + */ + doublereal m_sqrt_kbt; + + //! Current value of the pressure + doublereal m_press; + + //! Solution of the flux system + Array2D m_flux; + + //! saved value of the mixture thermal conductivity + doublereal m_lambda; + + //! Saved value of the mixture viscosity + doublereal m_viscmix; + + // work space + vector_fp m_spwork; + + //! Internal Function + + //! 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 updateViscosity_T(); + + //! Update the temperature-dependent parts of the mixture-averaged + //! thermal conductivity. + void updateCond_T(); + + //! Update the species viscosities + /*! + * Internal routine is run whenever the update_boolean + * m_spvisc_ok is false. This routine will calculate + * internal values for the species viscosities. + * + * @internal + */ + void updateSpeciesViscosities(); + + //! Update the binary diffusion coefficients wrt T. + /*! + * These are evaluated + * from the polynomial fits at unit pressure (1 Pa). + */ + void updateDiff_T(); + + //! Boolean indicating that mixture viscosity is current + bool m_viscmix_ok; + + //! Boolean indicating that weight factors wrt viscosity is current + bool m_viscwt_ok; + + //! Flag to indicate that the pure species viscosities + //! are current wrt the temperature + bool m_spvisc_ok; + + //! Boolean indicating that mixture diffusion coeffs are current + bool m_diffmix_ok; + + //! Boolean indicating that binary diffusion coeffs are current + bool m_bindiff_ok; + + //! Flag to indicate that the pure species conductivities + //! are current wrt the temperature + bool m_spcond_ok; + + //! Boolean indicating that mixture conductivity is current + bool m_condmix_ok; + + //! Mode for fitting the species viscosities + /*! + * Either its CK_Mode or its cantera mode + * in CK_Mode visc is fitted to a polynomial + * in Cantera mode sqrt(visc) is fitted. + */ + int m_mode; + + //! Internal storage for the diameter - diameter + //! species interactions + DenseMatrix m_diam; + + //! Debugging flags + /*! + * Turn on to get debugging information + */ + bool m_debug; + + //! Number of dimensions + /*! + * Either 1, 2, or 3 + */ + int m_nDim; + }; +} +#endif + + + + + + diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp new file mode 100644 index 000000000..07143106e --- /dev/null +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -0,0 +1,722 @@ +/** + * @file LiquidTransport.cpp + * Mixture-averaged transport properties for ideal gas mixtures. + */ +/* + * $Revision$ + * $Date$ + */ + +#include "ThermoPhase.h" +#include "LiquidTransport.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-20 + + +namespace Cantera { + + //////////////////// class LiquidTransport methods ////////////// + + + LiquidTransport::LiquidTransport() : + m_nsp(0), + m_tmin(-1.0), + m_tmax(100000.), + m_iStateMF(-1), + m_temp(-1.0), + m_logt(0.0), + m_sqrt_t(-1.0), + m_t14(-1.0), + m_t32(-1.0), + m_sqrt_kbt(-1.0), + m_press(-1.0), + m_lambda(-1.0), + m_viscmix(-1.0), + m_viscmix_ok(false), + m_viscwt_ok(false), + m_spvisc_ok(false), + m_diffmix_ok(false), + m_bindiff_ok(false), + m_spcond_ok(false), + m_condmix_ok(false), + m_mode(-1000), + m_debug(false), + m_nDim(1) + { + + + } + + + // Initialize the object + /* + * This is where we dimension everything. + */ + bool LiquidTransport::init(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()); + + // copy polynomials and parameters into local storage + m_poly = tr.poly; + viscCoeffsVector_ = tr.viscCoeffsVector_; + m_condcoeffs = tr.condcoeffs; + m_diffcoeffs = tr.diffcoeffs; + + m_mode = tr.mode; + m_diam = tr.diam; + m_eps = tr.eps; + m_alpha = tr.alpha; + + m_phi.resize(m_nsp, m_nsp, 0.0); + + + m_wratjk.resize(m_nsp, m_nsp, 0.0); + m_wratkj1.resize(m_nsp, m_nsp, 0.0); + int j, k; + for (j = 0; j < m_nsp; j++) + for (k = j; k < m_nsp; k++) { + m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]); + m_wratjk(k,j) = sqrt(m_wratjk(j,k)); + m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]); + } + + m_polytempvec.resize(5); + m_visc.resize(m_nsp); + m_sqvisc.resize(m_nsp); + m_cond.resize(m_nsp); + m_bdiff.resize(m_nsp, 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); + m_Grad_mu.resize(m_nDim * m_nsp, 0.0); + + + // set all flags to false + m_viscmix_ok = false; + m_viscwt_ok = false; + m_spvisc_ok = false; + m_spcond_ok = false; + m_condmix_ok = false; + m_spcond_ok = false; + m_diffmix_ok = false; + + return true; + } + + + + /****************** viscosity ******************************/ + + /* + * The viscosity is computed using the Wilke mixture rule. + * \f[ + * \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}. + * \f] + * Here \f$ \mu_k \f$ is the viscosity of pure species \e k, + * and + * \f[ + * \Phi_{k,j} = \frac{\left[1 + * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2} + * {\sqrt{8}\sqrt{1 + M_k/M_j}} + * \f] + * @see updateViscosity_T(); + */ + doublereal LiquidTransport::viscosity() { + + update_T(); + update_C(); + + if (m_viscmix_ok) return m_viscmix; + + // update m_visc[] and m_phi[] if necessary + if (!m_viscwt_ok) updateViscosity_T(); + + multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork)); + + m_viscmix = 0.0; + for (int k = 0; k < m_nsp; k++) { + m_viscmix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom; + } + return m_viscmix; + } + + + /******************* binary diffusion coefficients **************/ + + + void LiquidTransport::getBinaryDiffCoeffs(int ld, doublereal* d) { + int i,j; + + update_T(); + + // if necessary, evaluate the binary diffusion coefficents + // from the polynomial fits + if (!m_bindiff_ok) updateDiff_T(); + doublereal pres = m_thermo->pressure(); + + doublereal rp = 1.0/pres; + for (i = 0; i < m_nsp; i++) + for (j = 0; j < m_nsp; j++) { + d[ld*j + i] = rp * m_bdiff(i,j); + } + } + + + void LiquidTransport::getMobilities(doublereal* 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 LiquidTransport::set_Grad_V(const doublereal* grad_V) { + for (int a = 0; a < m_nDim; a++) { + m_Grad_V[a] = grad_V[a]; + } + } + + void LiquidTransport::set_Grad_T(const doublereal* grad_T) { + for (int a = 0; a < m_nDim; a++) { + m_Grad_T[a] = grad_T[a]; + } + } + + void LiquidTransport::set_Grad_X(const doublereal* grad_X) { + int itop = m_nDim * m_nsp; + for (int i = 0; i < itop; i++) { + m_Grad_X[i] = grad_X[i]; + } + } + + + /****************** thermal conductivity **********************/ + + /* + * The thermal conductivity is computed from the following mixture rule: + * \[ + * \lambda = 0.5 \left( \sum_k X_k \lambda_k + * + \frac{1}{\sum_k X_k/\lambda_k}\right) + * \] + */ + doublereal LiquidTransport::thermalConductivity() { + int k; + + update_T(); + update_C(); + + if (!m_spcond_ok) updateCond_T(); + if (!m_condmix_ok) { + doublereal sum1 = 0.0, sum2 = 0.0; + for (k = 0; k < m_nsp; k++) { + sum1 += m_molefracs[k] * m_cond[k]; + sum2 += m_molefracs[k] / m_cond[k]; + } + m_lambda = 0.5*(sum1 + 1.0/sum2); + } + return m_lambda; + } + + + /****************** thermal diffusion coefficients ************/ + + /** + * 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 LiquidTransport::getThermalDiffCoeffs(doublereal* dt) { + int k; + for (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 LiquidTransport::getSpeciesFluxes(int ndim, + 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); + } + + /** + * @param ndim The number of spatial dimensions (1, 2, or 3). + * @param grad_T The temperature gradient (ignored in this model). + * @param ldx Leading dimension of the grad_X array. + * The diffusive mass flux of species \e k is computed from + * + * \f[ + * \vec{j}_k = -n M_k D_k \nabla X_k. + * \f] + */ + void LiquidTransport::getSpeciesFluxesExt(int ldf, doublereal* fluxes) { + int n, k; + + 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]. + * + * For the single species case or the pure fluid case + * the routine returns the self-diffusion coefficient. + * This is need to avoid a Nan result in the formula + * below. + */ + void LiquidTransport::getMixDiffCoeffs(doublereal* d) { + + update_T(); + update_C(); + + // update the binary diffusion coefficients if necessary + if (!m_bindiff_ok) updateDiff_T(); + + int k, j; + doublereal mmw = m_thermo->meanMolecularWeight(); + doublereal sumxw = 0.0, sum2; + doublereal p = m_press; + if (m_nsp == 1) { + d[0] = m_bdiff(0,0) / p; + } else { + for (k = 0; k < m_nsp; k++) sumxw += m_molefracs[k] * m_mw[k]; + for (k = 0; k < m_nsp; k++) { + sum2 = 0.0; + for (j = 0; j < m_nsp; j++) { + if (j != k) { + sum2 += m_molefracs[j] / m_bdiff(j,k); + } + } + if (sum2 <= 0.0) { + d[k] = m_bdiff(k,k) / p; + } else { + d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2); + } + } + } + } + + + // 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 + */ + void LiquidTransport::update_T() + { + doublereal t = m_thermo->temperature(); + if (t == m_temp) return; + if (t < 0.0) { + throw CanteraError("LiquidTransport::update_T", + "negative temperature "+fp2str(t)); + } + + // Compute various functions of temperature + m_temp = t; + m_logt = log(m_temp); + m_kbt = Boltzmann * m_temp; + m_sqrt_t = sqrt(m_temp); + m_t14 = sqrt(m_sqrt_t); + m_t32 = m_temp * m_sqrt_t; + m_sqrt_kbt = sqrt(Boltzmann*m_temp); + + // compute powers of log(T) + m_polytempvec[0] = 1.0; + m_polytempvec[1] = m_logt; + m_polytempvec[2] = m_logt*m_logt; + m_polytempvec[3] = m_logt*m_logt*m_logt; + m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt; + + // temperature has changed, so polynomial temperature + // interpolations will need to be reevaluated. + // Set all of these flags to false + m_viscmix_ok = false; + m_spvisc_ok = false; + m_viscwt_ok = false; + m_spcond_ok = false; + m_diffmix_ok = false; + m_bindiff_ok = false; + m_condmix_ok = false; + + // For now, for a concentration redo also + m_iStateMF = -1; + } + + /** + * @internal This is called the first time any transport property + * is requested from Mixture after the concentrations + * have changed. + */ + void LiquidTransport::update_C() + { + + doublereal pres = m_thermo->pressure(); + // Check for changes in the mole fraction vector. + int iStateNew = m_thermo->getIStateMF(); + if (iStateNew == m_iStateMF) { + if (pres == m_press) { + return; + } + } else { + m_iStateMF = iStateNew; + } + m_press = pres; + + // signal that concentration-dependent quantities will need to + // be recomputed before use, and update the local mole + // fractions. + + m_viscmix_ok = false; + m_diffmix_ok = false; + m_condmix_ok = false; + + m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); + + // add an offset to avoid a pure species condition or + // negative mole fractions. MIN_X is 1.0E-20, a value + // which is below the additive machine precision of mole fractions. + int k; + for (k = 0; k < m_nsp; k++) { + m_molefracs[k] = fmaxx(MIN_X, m_molefracs[k]); + } + } + + + /************************************************************************* + * + * methods to update temperature-dependent properties + * + *************************************************************************/ + + /** + * Update the temperature-dependent parts of the mixture-averaged + * thermal conductivity. + */ + void LiquidTransport::updateCond_T() { + + int k; + if (m_mode == CK_Mode) { + for (k = 0; k < m_nsp; k++) { + m_cond[k] = exp(dot4(m_polytempvec, m_condcoeffs[k])); + } + } + else { + for (k = 0; k < m_nsp; k++) { + m_cond[k] = m_sqrt_t*dot5(m_polytempvec, m_condcoeffs[k]); + } + } + m_spcond_ok = true; + m_condmix_ok = false; + } + + + /** + * Update the binary diffusion coefficients. These are evaluated + * from the polynomial fits at unit pressure (1 Pa). + */ + void LiquidTransport::updateDiff_T() { + + // evaluate binary diffusion coefficients at unit pressure + int i,j; + int ic = 0; + if (m_mode == CK_Mode) { + for (i = 0; i < m_nsp; i++) { + for (j = i; j < m_nsp; j++) { + m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic])); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } + else { + for (i = 0; i < m_nsp; i++) { + for (j = i; j < m_nsp; j++) { + m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec, + m_diffcoeffs[ic]); + m_bdiff(j,i) = m_bdiff(i,j); + ic++; + } + } + } + + m_bindiff_ok = true; + m_diffmix_ok = false; + } + + + /** + * Update the pure-species viscosities. + */ + void LiquidTransport::updateSpeciesViscosities() { + + int k; + if (m_mode == CK_Mode) { + for (k = 0; k < m_nsp; k++) { + m_visc[k] = exp(dot4(m_polytempvec, viscCoeffsVector_[k])); + m_sqvisc[k] = sqrt(m_visc[k]); + } + } + else { + for (k = 0; k < m_nsp; k++) { + // the polynomial fit is done for sqrt(visc/sqrt(T)) + m_sqvisc[k] = m_t14*dot5(m_polytempvec, viscCoeffsVector_[k]); + m_visc[k] = (m_sqvisc[k]*m_sqvisc[k]); + } + } + m_spvisc_ok = true; + } + + + /** + * 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 LiquidTransport::updateViscosity_T() { + doublereal vratiokj, wratiojk, factor1; + + if (!m_spvisc_ok) updateSpeciesViscosities(); + + // see Eq. (9-5.15) of Reid, Prausnitz, and Poling + int j, k; + for (j = 0; j < m_nsp; j++) { + for (k = j; k < m_nsp; k++) { + vratiokj = m_visc[k]/m_visc[j]; + wratiojk = m_mw[j]/m_mw[k]; + + // Note that m_wratjk(k,j) holds the square root of + // m_wratjk(j,k)! + factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j); + m_phi(k,j) = factor1*factor1 / + (SqrtEight * m_wratkj1(j,k)); + m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk); + } + } + m_viscwt_ok = true; + } + + /** + * This function returns a Transport data object for a given species. + * + */ + struct GasTransportData LiquidTransport:: + getGasTransportData(int kSpecies) + { + struct GasTransportData td; + td.speciesName = m_thermo->speciesName(kSpecies); + + + td.wellDepth = m_eps[kSpecies] / Boltzmann; + td.diameter = m_diam(kSpecies, kSpecies) * 1.0E10; + td.polarizability = m_alpha[kSpecies] * 1.0E30; + + + return td; + } + + + /* + * + * Solve for the diffusional velocities in the Stefan-Maxwell equations + * + */ + void LiquidTransport::stefan_maxwell_solve() { + int i, j, a; + + int VIM = 2; + m_B.resize(m_nsp, VIM); + //! grab a local copy of the molecular weights + const vector_fp& M = m_thermo->molecularWeights(); + + + //! get the mean molecular weight of the mixture + //double M_mix = m_thermo->meanMolecularWeight(); + + + //! get the concentration of the mixture + //double rho = m_thermo->density(); + //double c = rho/M_mix; + + + m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); + + double T = m_thermo->temperature(); + + + /* electrochemical potential gradient */ + for (i = 0; i < m_nsp; i++) { + for (a = 0; a < VIM; a++) { + m_Grad_mu[a*m_nsp + i] = m_chargeSpecies[i] * Faraday * m_Grad_V[a] + + (GasConstant*T/m_molefracs[i]) * m_Grad_X[a*m_nsp+i]; + } + } + + /* + * Just for Note, m_A(i,j) refers to the ith row and jth column. + * They are still fortran ordered, so that i varies fastest. + */ + switch ( VIM ) { + case 1: /* 1-D approximation */ + m_B(0,0) = 0.0; + for (j = 0; j < m_nsp; j++) { + m_A(0,j) = 1.0; + } + for (i = 1; i < m_nsp; i++){ + m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T); + for (j = 0; j < m_nsp; j++){ + if (j != i) { + m_A(i,j) = m_molefracs[i] / ( M[j] * m_DiffCoeff_StefMax(i,j)); + m_A(i,i) -= m_molefracs[j] / ( M[i] * m_DiffCoeff_StefMax(i,j)); + } + else if (j == i) { + m_A(i,i) = 0.0; + } + } + } + + //! invert and solve the system Ax = b. Answer is in m_B + solve(m_A, m_B.ptrColumn(0)); + + m_flux = m_B; + + + break; + case 2: /* 2-D approximation */ + m_B(0,0) = 0.0; + m_B(0,1) = 0.0; + for (j = 0; j < m_nsp; j++) { + m_A(0,j) = 1.0; + } + for (i = 1; i < m_nsp; i++){ + m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T); + m_B(i,1) = m_concentrations[i] * m_Grad_mu[m_nsp + i] / (GasConstant * T); + for (j = 0; j < m_nsp; j++){ + if (j != i) { + m_A(i,j) = m_molefracs[i] / ( M[j] * m_DiffCoeff_StefMax(i,j)); + m_A(i,i) -= m_molefracs[j] / ( M[i] * m_DiffCoeff_StefMax(i,j)); + } + else if (j == i) { + m_A(i,i) = 0.0; + } + } + } + + //! invert and solve the system Ax = b. Answer is in m_B + solve(m_A, m_B); + + m_flux = m_B; + + + break; + + case 3: /* 3-D approximation */ + m_B(0,0) = 0.0; + m_B(0,1) = 0.0; + m_B(0,2) = 0.0; + for (j = 0; j < m_nsp; j++) { + m_A(0,j) = 1.0; + } + for (i = 1; i < m_nsp; i++){ + m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T); + m_B(i,1) = m_concentrations[i] * m_Grad_mu[m_nsp + i] / (GasConstant * T); + m_B(i,2) = m_concentrations[i] * m_Grad_mu[2*m_nsp + i] / (GasConstant * T); + for (j = 0; j < m_nsp; j++){ + if (j != i) { + m_A(i,j) = m_molefracs[i] / ( M[j] * m_DiffCoeff_StefMax(i,j)); + m_A(i,i) -= m_molefracs[j] / ( M[i] * m_DiffCoeff_StefMax(i,j)); + } + else if (j == i) { + m_A(i,i) = 0.0; + } + } + } + + //! invert and solve the system Ax = b. Answer is in m_B + solve(m_A, m_B); + + m_flux = m_B; + + + break; + default: + printf("uninmplemetnd\n"); + throw CanteraError("routine", "not done"); + break; + } + + + } +} diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h new file mode 100644 index 000000000..695088ea8 --- /dev/null +++ b/Cantera/src/transport/LiquidTransport.h @@ -0,0 +1,634 @@ +/** + * + * @file LiquidTransport.h + * Header file defining class LiquidTransport + */ +/* + * $Revision$ + * $Date$ + */ + +#ifndef CT_LIQUIDTRAN_H +#define CT_LIQUIDTRAN_H + + + +// STL includes +#include +#include +#include +#include +#include + +using namespace std; + +// Cantera includes +#include "TransportBase.h" +#include "DenseMatrix.h" +#include "LiquidTransportParams.h" + +namespace Cantera { + + const int LVISC_CONSTANT = 0; + const int LVISC_WILKES = 1; + + + + class TransportParams; + + + //! 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 LiquidTransport : public Transport { + + public: + + //! virtual destructor + virtual ~LiquidTransport() {} + + //! Return the model id for this transport parameterization + virtual int model() { return cLiquidTransport; } + + //! overloaded base class methods + + //! Returns the viscosity of the solution + /*! + * The viscosity is computed using the Wilke mixture rule. + * \f[ + * \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}. + * \f] + * Here \f$ \mu_k \f$ is the viscosity of pure species \e k, + * and + * \f[ + * \Phi_{k,j} = \frac{\left[1 + * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2} + * {\sqrt{8}\sqrt{1 + M_k/M_j}} + * \f] + * @see updateViscosity_T(); + * + * Controlling update boolean m_viscmix_ok + */ + virtual doublereal viscosity(); + + //! Returns the pure species viscosities + /*! + * + * Controlling update boolean = m_viscwt_ok + */ + virtual void getSpeciesViscosities(doublereal* visc) + { updateViscosity_T(); copy(m_visc.begin(), m_visc.end(), visc); } + + virtual void getThermalDiffCoeffs(doublereal* dt); + + //! Return the thermal conductivity of the solution + /*! + * The thermal conductivity is computed from the following mixture rule: + * \f[ + * \lambda = 0.5 \left( \sum_k X_k \lambda_k + * + \frac{1}{\sum_k X_k/\lambda_k}\right) + * \f] + * + * Controlling update boolean = m_condmix_ok + */ + virtual doublereal thermalConductivity(); + + //! Returns the binary diffusion coefficients + /*! + * @param ld + * @param d + */ + virtual void getBinaryDiffCoeffs(int ld, doublereal* 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* d); + + + //! Get the Mobilities + /*! + * @param mobil + */ + virtual void getMobilities(doublereal* mobil); + + //! Specify the value of the gradient of the voltage + /*! + * + * @param grad_V Gradient of the voltage (length num dimensions); + */ + virtual void set_Grad_V(const doublereal* 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* 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* grad_X); + + //! 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 + */ + virtual void update_T(); + + //! Handles the effects of changes in the mixture concentration + /*! + * This is called the first time any transport property + * is requested from Mixture after the concentrations + * have changed. + * + * @internal + */ + virtual void update_C(); + + /** + * @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, + doublereal* grad_T, + int ldx, const doublereal* grad_X, + int ldf, doublereal* 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 getSpeciesFluxesExt(int ldf, doublereal* fluxes); + + + //! 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 init(LiquidTransportParams& tr); + + friend class TransportFactory; + + /** + * Return a structure containing all of the pertinent parameters + * about a species that was used to construct the Transport + * properties in this object. + * + * @param k Species number to obtain the properties about. + */ + struct GasTransportData getGasTransportData(int k); + + + //! Solve the stefan_maxell equations for the diffusive fluxes. + void stefan_maxwell_solve(); + protected: + + //! default constructor + LiquidTransport(); + + private: + + + + //! Number of species in the mixture + int m_nsp; + + //! 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; + + // polynomial fits + vector > m_poly; + + //! Polynomial coefficients of the viscosity + /*! + * These express the temperature dependendence of the pures + * species viscosities. + */ + std::vector viscCoeffsVector_; + + //! Polynomial coefficients of the conductivities + /*! + * These express the temperature dependendence of the pures + * species conductivities + */ + vector m_condcoeffs; + + //! Polynomial coefficients of the binary diffusion coefficients + /*! + * These express the temperature dependendence of the + * binary diffusivities. An overall pressure dependence is then + * added. + */ + vector m_diffcoeffs; + + + //! Internal value of the gradient of the mole fraction vector + /*! + * 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 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; + + //! Gradient of the electrochemical potential + /*! + * m_nsp is the number of species in the fluid + * k is the species index + * n is the dimensional index (x, y, or z) + * + * m_Grad_mu[n*m_nsp + k] + */ + vector_fp m_Grad_mu; + + // property values + + //! Array of Binary Diffusivities + /*! + * This has a size equal to nsp x nsp + * It is a symmetric matrix. + * D_ii is undefined. + * + * units m2/sec + */ + DenseMatrix m_bdiff; + + //! Species viscosities + /*! + * Viscosity of the species + * Length = number of species + * + * Depends on the temperature and perhaps pressure, but + * not the species concentrations + * + * controlling update boolean -> m_spvisc_ok + */ + vector_fp m_visc; + + //! Sqrt of the species viscosities + /*! + * The sqrt(visc) is used in the mixing formulas + * Length = m_nsp + * + * Depends on the temperature and perhaps pressure, but + * not the species concentrations + * + * controlling update boolean m_spvisc_ok + */ + vector_fp m_sqvisc; + + //! 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_spcond_ok + */ + vector_fp m_cond; + + //! Polynomials of the log of the temperature + vector_fp m_polytempvec; + + //! State of the mole fraction vector. + int m_iStateMF; + + //! Local copy of the mole fractions of the species in the phase + /*! + * Update info? + * length = m_nsp + */ + vector_fp m_molefracs; + + //! Local copy of the concentrations of the species in the phase + /*! + * Update info? + * length = m_nsp + */ + vector_fp m_concentrations; + + //! Local copy of the charge of each species + /*! + * Contains the charge of each species (length m_nsp) + */ + vector_fp m_chargeSpecies; + + //! Stefan-Maxwell Diffusion Coefficients at T, P and C + /*! + * These diffusion coefficients are considered to be + * a function of Temperature, Pressure, and Concentration. + */ + DenseMatrix m_DiffCoeff_StefMax; + + //! Viscosity model + /*! + * + */ + int viscosityModel_; + + //! viscosity weighting functions + DenseMatrix m_phi; + + //! Matrix of the ratios of the species molecular weights + /*! + * m_wratjk(i,j) = (m_mw[j]/m_mw[k])**0.25 + */ + DenseMatrix m_wratjk; + + //! Matrix of the ratios of the species molecular weights + /*! + * m_wratkj1(i,j) = (1.0 + m_mw[k]/m_mw[j])**0.5 + */ + DenseMatrix m_wratkj1; + + //! RHS to the stefan-maxwell equation + Array2D m_B; + + //! Matrix for the stefan maxwell equation. + DenseMatrix m_A; + + //! Internal storage for the species LJ well depth + vector_fp m_eps; + + //! Internal storage for species polarizability + vector_fp m_alpha; + + + + //! Current Temperature -> locally storred + /*! + * This is used to test whether new temperature computations + * should be performed. + */ + doublereal m_temp; + + //! Current log(T) + doublereal m_logt; + + //! Current value of kT + doublereal m_kbt; + + //! Current Temperature **0.5 + doublereal m_sqrt_t; + + //! Current Temperature **0.25 + doublereal m_t14; + + //! Current Temperature **1.5 + doublereal m_t32; + + //! Current temperature function + /*! + * This is equal to sqrt(Boltzmann * T) + */ + doublereal m_sqrt_kbt; + + //! Current value of the pressure + doublereal m_press; + + //! Solution of the flux system + Array2D m_flux; + + //! saved value of the mixture thermal conductivity + doublereal m_lambda; + + //! Saved value of the mixture viscosity + doublereal m_viscmix; + + // work space + vector_fp m_spwork; + + //! Internal Function + + //! 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 updateViscosity_T(); + + //! Update the temperature-dependent parts of the mixture-averaged + //! thermal conductivity. + void updateCond_T(); + + //! Update the species viscosities + /*! + * Internal routine is run whenever the update_boolean + * m_spvisc_ok is false. This routine will calculate + * internal values for the species viscosities. + * + * @internal + */ + void updateSpeciesViscosities(); + + //! Update the binary diffusion coefficients wrt T. + /*! + * These are evaluated + * from the polynomial fits at unit pressure (1 Pa). + */ + void updateDiff_T(); + + //! Boolean indicating that mixture viscosity is current + bool m_viscmix_ok; + + //! Boolean indicating that weight factors wrt viscosity is current + bool m_viscwt_ok; + + //! Flag to indicate that the pure species viscosities + //! are current wrt the temperature + bool m_spvisc_ok; + + //! Boolean indicating that mixture diffusion coeffs are current + bool m_diffmix_ok; + + //! Boolean indicating that binary diffusion coeffs are current + bool m_bindiff_ok; + + //! Flag to indicate that the pure species conductivities + //! are current wrt the temperature + bool m_spcond_ok; + + //! Boolean indicating that mixture conductivity is current + bool m_condmix_ok; + + //! Mode for fitting the species viscosities + /*! + * Either its CK_Mode or its cantera mode + * in CK_Mode visc is fitted to a polynomial + * in Cantera mode sqrt(visc) is fitted. + */ + int m_mode; + + //! Internal storage for the diameter - diameter + //! species interactions + DenseMatrix m_diam; + + //! Debugging flags + /*! + * Turn on to get debugging information + */ + bool m_debug; + + //! Number of dimensions + /*! + * Either 1, 2, or 3 + */ + int m_nDim; + }; +} +#endif + + + + + + diff --git a/Cantera/src/transport/LiquidTransportParams.h b/Cantera/src/transport/LiquidTransportParams.h new file mode 100644 index 000000000..b2afd4368 --- /dev/null +++ b/Cantera/src/transport/LiquidTransportParams.h @@ -0,0 +1,62 @@ +#ifndef CT_LIQUIDTRANSPORTPARAMS_H +#define CT_LIQUIDTRANSPORTPARAMS_H + +#include + +#include "ct_defs.h" +#include "TransportBase.h" +#include "xml.h" +#include "XML_Writer.h" + +namespace Cantera { + + /** + * + * Holds transport data. Used by TransportFactory. + * + */ + class LiquidTransportParams { + + public: + + LiquidTransportParams() : thermo(0), xml(0) {} + virtual ~LiquidTransportParams(); + int nsp; + + // phase_t* mix; + thermo_t* thermo; + vector_fp mw; + + // polynomial fits + std::vector viscCoeffsVector_; + std::vector condcoeffs; + std::vector diffcoeffs; + vector_fp polytempvec; + + std::vector > poly; + std::vector omega22_poly; + std::vector astar_poly; + std::vector bstar_poly; + std::vector cstar_poly; + + vector_fp zrot; + vector_fp crot; + + std::vector polar; + vector_fp alpha; + vector_fp fitlist; + vector_fp eps; + vector_fp sigma; + DenseMatrix reducedMass; + DenseMatrix diam; + DenseMatrix epsilon; + DenseMatrix dipole; + DenseMatrix delta; + doublereal tmax, tmin; + int mode; + XML_Writer* xml; + int log_level; + }; +} + +#endif diff --git a/Cantera/src/transport/Makefile.in b/Cantera/src/transport/Makefile.in index d61a6dfe4..caf985c79 100644 --- a/Cantera/src/transport/Makefile.in +++ b/Cantera/src/transport/Makefile.in @@ -32,10 +32,11 @@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG) # Transport Object Files OBJS = TransportFactory.o MultiTransport.o MixTransport.o MMCollisionInt.o \ - SolidTransport.o DustyGasTransport.o + SolidTransport.o DustyGasTransport.o LiquidTransport.o + TRAN_H = TransportFactory.h MultiTransport.h MixTransport.h \ MMCollisionInt.h SolidTransport.h DustyGasTransport.h \ - TransportBase.h L_matrix.h TransportParams.h + TransportBase.h L_matrix.h TransportParams.h LiquidTransport.h CXX_INCLUDES = -I../base -I../thermo -I../numerics @CXX_INCLUDES@ LIB = @buildlib@/libtransport.a diff --git a/Cantera/src/transport/MixTransport.h b/Cantera/src/transport/MixTransport.h index 2640b3932..6b63db49f 100755 --- a/Cantera/src/transport/MixTransport.h +++ b/Cantera/src/transport/MixTransport.h @@ -112,6 +112,14 @@ namespace Cantera { const doublereal* grad_X, int ldf, doublereal* fluxes); + //! 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 init(TransportParams& tr); friend class TransportFactory; diff --git a/Cantera/src/transport/TransportBase.h b/Cantera/src/transport/TransportBase.h index 7c3ed7f25..dcab7ad54 100755 --- a/Cantera/src/transport/TransportBase.h +++ b/Cantera/src/transport/TransportBase.h @@ -43,6 +43,7 @@ namespace Cantera { const int cUserTransport = 500; const int cFtnTransport = 600; const int cLiquidTransport = 700; + const int cAqueousTransport = 750; const int cRadiativeTransport = 800; // forward reference diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index 342875e5b..b31ee702c 100755 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -21,6 +21,8 @@ #include "SolidTransport.h" #include "DustyGasTransport.h" +#include "LiquidTransport.h" +#include "AqueousTransport.h" #include "TransportFactory.h" @@ -299,6 +301,7 @@ namespace Cantera { if (transportModel == "") return new Transport; + vector_fp state; Transport *tr = 0, *gastr = 0; DustyGasTransport* dtr = 0; @@ -334,15 +337,24 @@ namespace Cantera { dtr = (DustyGasTransport*)tr; dtr->initialize(phase, gastr); break; + case cLiquidTransport: + tr = new LiquidTransport; + tr->setThermo(*phase); + break; + case cAqueousTransport: + tr = new AqueousTransport; + tr->setThermo(*phase); + break; default: throw CanteraError("newTransport","unknown transport model"); } phase->restoreState(state); return tr; } + - /** +/** * Prepare to build a new kinetic-theory-based transport manager * for low-density gases. Uses polynomial fits to Monchick & Mason * collision integrals.