Added stubs for liquid phase and aqueous phase transport.
This commit is contained in:
Harry Moffat
parent
8737bda5f9
commit
0127285b88
9 changed files with 2799 additions and 3 deletions
735
Cantera/src/transport/AqueousTransport.cpp
Normal file
735
Cantera/src/transport/AqueousTransport.cpp
Normal file
|
|
@ -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 <iostream>
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
621
Cantera/src/transport/AqueousTransport.h
Normal file
621
Cantera/src/transport/AqueousTransport.h
Normal file
|
|
@ -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 <vector>
|
||||
#include <string>
|
||||
#include <map>
|
||||
#include <numeric>
|
||||
#include <algorithm>
|
||||
|
||||
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.
|
||||
*
|
||||
* <H2> Viscosity Calculation </H2>
|
||||
*
|
||||
* 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<vector<int> > m_poly;
|
||||
|
||||
//! Polynomial coefficients of the viscosity
|
||||
/*!
|
||||
* These express the temperature dependendence of the pures
|
||||
* species viscosities.
|
||||
*/
|
||||
vector<vector_fp> m_visccoeffs;
|
||||
|
||||
//! Polynomial coefficients of the conductivities
|
||||
/*!
|
||||
* These express the temperature dependendence of the pures
|
||||
* species conductivities
|
||||
*/
|
||||
vector<vector_fp> 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<vector_fp> 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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
722
Cantera/src/transport/LiquidTransport.cpp
Normal file
722
Cantera/src/transport/LiquidTransport.cpp
Normal file
|
|
@ -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 <iostream>
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
634
Cantera/src/transport/LiquidTransport.h
Normal file
634
Cantera/src/transport/LiquidTransport.h
Normal file
|
|
@ -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 <vector>
|
||||
#include <string>
|
||||
#include <map>
|
||||
#include <numeric>
|
||||
#include <algorithm>
|
||||
|
||||
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.
|
||||
*
|
||||
* <H2> Viscosity Calculation </H2>
|
||||
*
|
||||
* 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<vector<int> > m_poly;
|
||||
|
||||
//! Polynomial coefficients of the viscosity
|
||||
/*!
|
||||
* These express the temperature dependendence of the pures
|
||||
* species viscosities.
|
||||
*/
|
||||
std::vector<vector_fp> viscCoeffsVector_;
|
||||
|
||||
//! Polynomial coefficients of the conductivities
|
||||
/*!
|
||||
* These express the temperature dependendence of the pures
|
||||
* species conductivities
|
||||
*/
|
||||
vector<vector_fp> 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<vector_fp> 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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
62
Cantera/src/transport/LiquidTransportParams.h
Normal file
62
Cantera/src/transport/LiquidTransportParams.h
Normal file
|
|
@ -0,0 +1,62 @@
|
|||
#ifndef CT_LIQUIDTRANSPORTPARAMS_H
|
||||
#define CT_LIQUIDTRANSPORTPARAMS_H
|
||||
|
||||
#include <vector>
|
||||
|
||||
#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<vector_fp> viscCoeffsVector_;
|
||||
std::vector<vector_fp> condcoeffs;
|
||||
std::vector<vector_fp> diffcoeffs;
|
||||
vector_fp polytempvec;
|
||||
|
||||
std::vector<std::vector<int> > poly;
|
||||
std::vector<vector_fp > omega22_poly;
|
||||
std::vector<vector_fp > astar_poly;
|
||||
std::vector<vector_fp > bstar_poly;
|
||||
std::vector<vector_fp > cstar_poly;
|
||||
|
||||
vector_fp zrot;
|
||||
vector_fp crot;
|
||||
|
||||
std::vector<bool> 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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue